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ABSTRACT 


Sea  ice  is  of  fundamental  importance  in  weather,  climate  and  other 
geophysical  processes.  It  is  also  an  important  factor  for  naval 
operations  in  the  polar  regions,  in  particular  regarding  transport  of 
personnel  and  material  in  regions  where  sea  ice  is  likely  to  be  found 
and  assessment  and  prediction  of  acoustic  environments  in  polar 
regions.  Because  sea  ice  has  a  large  geographic  extent  and  short  time 
scale  for  variability  synthetic  aperture  radar  (SAR)  is  a  valuable 
technique  in  studying  sea  ice,  particularly  since  images  can  be  collected 
through  clouds  and  at  night.  SAR  information  on  sea  ice  will  be 
available  from  several  satellites  in  the  1990's  (ERS-1,  JERS-1,  Almaz, 
Radarsat  and  possibly  an  EOS  SAR).  Automated  interpretation 
techniques  are  required  because  of  the  large  number  and  high 
information  content  of  the  SAR  images  becoming  available.  Here  we 
report  research  on  automated-computer-based  techniques  for  such 
interpretation.  The  general  problem  that  we  face  is  to  extract 
geophysical  information  from  one  or  more  SAR  images.  The  work 
reported  here  concerns  estimation  of  sea  ice  motion.  Hence,  we  collect 
two  images  of  a  given  geographical  area  some  several  days  apart  and 
note  the  change  in  sea  ice  location  which  gives  us  the  average  sea  ice 
velocity  vector.  While  one  can  do  this  manually  by  establishing  tie 
points  on  the  same  piece  of  ice  in  the  two  scenes,  an  automated 
algorithm  is  preferable  for  both  cost  and  consistency  reasons.  The 
research  reported  here  covers  the  image  pyramid  area  correlation 
(IPAC)  and  feature  tracking  algorithms.  These  two  algorithms,  when 
used  together,  can  provide  robust  ice  velocity  estimates  for  ice  that 
translates  and  rotates.  A  particularly  useful  result  is  an  error  checking 
algorithm  that  graphically  displays  the  regions  where  most  errors  occur. 
The  SAR  interpretation  algorithms  developed  under  this  research 
funding  have  contributed  to  the  Geophysical  Data  Processing  System 
(GPS)  at  the  Alaskan  SAR  Facility  (ASF)  in  Fairbanks,  Alaska.  The 
geophysical  data  flowing  from  this  system  are  in  use  in  the  Joint  Ice 
Center  (JIC)  as  well  as  the  Naval  Polar  Oceanography  Center  (NPOC)  and 
the  remote  sensing  activities  of  the  Naval  Research  Laboratory, 
Detachment  South.  The  SAR  interpretation  algorithms  to  which  this 
research  contributes  assure  that  SAR  image  information  is  available  in  a 
timely  manner  for  use  in  ice  science  and  naval  applications. 
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I.  Introduction  and  Motivation  for  Sea  Ice  Research 

Sea  ice  is  of  fundamental  importance  in  weather,  climate  and  other 
geophysical  processes.  It  is  also  an  important  factor  for  naval 
operations  in  the  polar  regions,  in  particular  regarding  transport  of 
personnel  and  material  in  regions  where  sea  ice  is  likely  to  be  found 
and  assessment  and  prediction  of  acoustic  environments  in  polar 
regions.  Because  sea  ice  has  a  large  geographic  extent  and  short  time 
scale  for  variability  synthetic  aperture  radar  (SAR)  is  a  valuable 
technique  in  studying  sea  ice.  A  SAR  image  is  a  'radar  picture'  of  a 
scene  on  the  Earth's  surface.  The  picture  elements  or  pixels  of  the  scene 
are  a  map  of  the  radar  reflectivity  of  the  pixel's  area  on  the  surface. 

SAR  images  are  particularly  useful  in  that  they  can  be  collected  through 
clouds  and  at  night.  Such  a  capability  is  clearly  important  in  polar 
regions.  SAR  information  on  sea  ice  will  be  available  from  several 
satellites  in  the  1990's  (ERS-1,  JERS-1,  Almaz,  Radarsat  and  possibly  an 
EOS  SAR).  Automated  interpretation  techniques  are  required  because  of 
the  large  number  and  high  information  content  of  the  SAR  images 
becoming  available.  We  have  worked  to  develop  automated-computer- 
based  techniques  for  such  interpretation  and  make  use  of  advanced 
methods  and  concepts  in  image  processing,  computer  vision  and 
artificial  intelligence. 

The  general  problem  that  we  face  is  to  extract  geophysical 
information  from  one  or  more  SAR  images.  The  work  reported  here 
concerns  estimation  of  sea  ice  motion.  Hence,  we  collect  two  images  of  a 
given  geographical  area  some  several  days  apart  and  note  the  change  in 
sea  ice  location  which  gives  us  the  average  sea  ice  velocity  vector. 

While  one  can  do  this  manually  by  establishing  tie  points  on  the  same 
piece  of  ice  in  the  two  scenes,  an  automated  algorithm  is  preferable  for 
both  cost  and  consistency  reasons.  The  research  reported  here  contains 
important  results  for  the  solution  of  this  problem. 

The  SAR  interpretation  algorithms  developed  under  this  research 
funding  have  have  contributed  to  the  Geophysical  Data  Processing 
System  (GPS)  at  the  Alaskan  SAR  Facility  (ASF)  in  Fairbanks,  Alaska. 

The  geophysical  data  flowing  from  this  system  are  in  use  in  the  Joint  Ice 
Center  (JIC)  as  well  as  the  Naval  Polar  Oceanography  Center  (NPOC)  and 
the  remote  sensing  activities  of  the  Naval  Research  Laboratory, 
Detachment  South.  The  SAR  interpretation  algorithms  to  which  this 
research  contributes  assure  that  SAR  image  information  is  available  in  a 
timely  manner  for  use  in  ice  science  and  naval  applications. 
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An  example  of  the  contribution  that  our  results  have  made  is 
illustrated  in  Fig.  6  on  p.  63  of  Vesecky  et  al.  (1988)  --  see  the 
Appendix,  beginning  on  page  12,  for  a  copy  of  this  paper.  Our  general 
scheme  for  ice  tracking,  illustrated  in  Fig.  6,  was  first  presented  in  1988 
and  eventually  lead  to  the  basis  for  the  algorithm  implemented  in  the 
Geophysical  Processor  at  the  ASF.  Ideas  in  the  results  reported  here 
were  fundamental  to  both  the  general  scheme  of  the  ice  tracking 
algorithm  as  well  as  the  components  within  the  algorithm. 


II.  Research  Objectives 

The  objectives  of  research  reported  here  follow  from  the 
circumstances  discussed  above.  The  objectives  can  be  summarized  as 
follows: 

1.  Investigate  methods  for  estimating  sea  ice  motion  given  two  more 
SAR  images  of  sea  ice  fields,  collected  at  intervals  of  the  order  of  days. 

2.  Develop  the  methods  above  into  working  computational  algorithms 
that  function  in  an  automated  fashion. 

3.  Test  the  algorithms  above  on  image  pairs  collected  by  Seasat  SAR 
and  a  variety  of  aircraft  SAR  images. 

5.  Modify  the  algorithms  of  item  2  in  the  light  of  3  until  satisfactory 
performance  is  achieved. 

6.  Support  the  implementation  of  the  algorithms  developed  on  the 
Geophysical  Processor  System  of  the  Alaskan  SAR  Facility. 


III.  Research  Results 

The  research  results  stemming  from  the  objectives  above  can  be 
summarized  under  four  topics: 

A.  Image  pyramid  area  correlation  schemes 

B.  Feature  tracing  schemes 

C.  Rotation  invariant  correlation 

D.  Checking  results  by  image  subtraction. 
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The  most  significant  progress  made  on  these  topics  concerns  the  first 
two.  The  results  in  these  two  areas  are  incorporated  in  the  sea  ice 
tracking  algorithms  now  in  place  at  the  Alaskan  SAR  Facility  at  the 
University  of  Alaska  in  Fairbanks. 

A.  Image  Pyramid  Area  Correlation  Schemes 

For  tracking  of  sea  ice  that  translates  without  significant  rotation 
the  image  pyramid  area  correlation  scheme  has  proved  to  be  a  robust 
method  of  ice  tracking.  In  this  scheme  a  sequence  of  images  or  image 
pyramid  is  constructed  in  which  each  image  in  the  sequence  has 
progressively  cruder  resolution.  A  typical  scheme  is  to  average  4  pixels 
in  one  scene  to  produce  the  next  image  in  the  sequence.  The  pyramid 
refers  to  the  fact  that  as  resolution  is  degraded  the  images  have  fewer 
pixels  (covering  more  surface  area  each)  and  the  images  become  smaller 
in  terms  of  pixels  although  the  image  corresponds  to  the  same 
geographical  area.  These  progressively  smaller  images  can  be  stacked 
up  to  form  a  pyramid  as  shown  in  Fig.  1  of  Vesecky  et  al.  ( 1988a)  in  the 
Appendix. 

The  image  pyramid  area  correlation  method  for  ice  tracking  begins 
with  the  crude  resolution  image  pair  at  the  top  of  the  pyramid.  A  small 
square  block  of  pixels,  a  subimage,  in  the  first  image,  image  1,  of  a  pair 
is  selected  and  cross  correlated  with  areas  in  the  second  image  of  the 
pair,  image  2.  The  region  with  highest  correlation  in  image  2  shows 
where  the  ice  in  the  selected  subimage  has  moved  over  the  several  days 
between  the  collection  of  image  1  and  image  2.  This  process  is  repeated 
for  all  the  non-overlapping  subimages  of  image  1.  The  result  is  a  field 
of  vectors  showing  a  crude  estimate  of  the  large  scale  ice  motion  field. 
The  process  is  repeated  with  the  next  image  pair  in  the  pyramid,  i.e. 
with  slightly  better  resolution.  However,  we  now  have  a  general  idea  of 
the  large  scale  motion  field;  so  we  know  approximately  where  to 
conduct  our  correlation  search  -  we  don't  have  to  cross  correlate  with 
all  of  image  2,  just  a  suitable,  small  part  of  it.  After  the  process  has 
been  repeated  several  times  we  arrive  at  the  base  of  the  image 
pyramid,  obtaining  the  highest  resolution  velocity  vector  field.  At  each 
step  we  have  used  the  velocity  field  estimate  at  the  next  higher  level  to 
guide  our  search  for  the  highest  cross  correlation  match.  Results  of  this 
algorithm  are  shown  in  Fig.  2  of  Vesecky  et  al.  ( 1988a)  or  Fig.  4.  of 
Vesecky  et  al.  (1988b). 
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B.  Feature  Tracking  Schemes 

While  the  image  pyramid  area  correlation  method  works  very 
robustly  on  ice  that  translates,  but  does  not  rotate,  it  begins  to  break 
down  when  the  ice  in  the  image  rotates  or  shears.  As  shown  in  Vesecky 
et  al.  (1988a),  a  rotation  of  only  7°  reduces  correlation  in  a  typical  SAR 
sea  ice  image  to  0.8.  A  shear  of  about  10°  is  required  to  reduce  the 
correlation  to  0.8.  Thus,  even  a  small  amount  of  rotation  or  shear  can 
cause  mistakes  in  the  image  pyramid  method.  Regions  where  such 
errors  are  likely  to  occur  can  be  recognized  by  the  fact  that  the 
maximum  cross  correlation  is  reduced  below  0.8.  In  such  regions  it  is 
necessary  to  use  another  method  or  to  expand  the  IPAC  scheme  to 
account  for  rotation.  Expanding  the  IPAC  to  include  all  possible 
rotations  would  mean  expanding  the  search  space  by  a  factor  of  50  or 
more  with  the  associated  increase  in  computation  load.  Hence,  a  second 
method  appears  attractive. 

Feature  tracking  seeks  to  decompose  each  SAR  sea  ice  image  into  a 
collection  of  features.  Vesecky  et  al.  ( 1988b)  did  this  for  the  first  time. 
The  features  used  were  segments  of  the  edges  of  floes  as  shown  in  Fig.  9 
of  the  above  reference  (reproduced  in  the  Appendix).  The  trick  here  is 
that  the  features  are  characterized  independent  of  their  orientation. 

One  then  seeks  to  find  a  matching  pair  of  features  in  both  images.  Once 
a  matching  pair  of  features  is  found  both  the  translation  and  rotation  of 
the  ice  are  known.  While  the  feature  tracking  scheme  does  not  yield  a 
complete  ice  velocity  map,  it  does  show  portions  of  the  map  indicating 
where  rotation  occurs.  Thus,  by  using  feature  tracking  one  can 
determine  the  amount  of  translation  and  rotation  in  a  particular  region 
and  adapt  IPAC  by  knowing  the  amount  of  rotation  to  apply  when 
searching  for  corresponding  regions  in  an  image  pair. 

A  useful  approach  to  tracking  sea  ice  that  rotates  is  to  use  a 
combination  of  IPAC  and  feature  tracking,  as  suggested  in  Vesecky  et  al. 
(1988a).  For  example,  one  could  first  run  the  feature  tracking 
algorithm,  noting  regions  of  significant  rotation.  This  information  would 
be  used  as  part  of  the  guided  search  for  matching  regions  in  IPAC.  An 
alternative  would  be  to  run  IPAC  first  and  then  seek  feature  tracking 
information  only  in  those  regions  where  no  high  correlation  matches 
were  found.  In  these  'no  match'  cases  feature  tracking  would  be  done 
to  provide  rotation  information  so  that  a  match  could  be  found. 

Schemes  such  as  these  are  now  used  in  the  geophysical  processor  at  the 
Alaskan  SAR  Facility  (ASF)  in  Fairbanks. 
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C.  Rotation  Invarient  Correlation 

Another  approach  to  tracking  ice  that  rotates  is  rotation  invariant 
correlation.  Vesecky  et  al.  (1987)  applied  the  rotation  invariant 
moments  of  Hu  (1972)  to  the  tracking  of  rotating  sea  ice  floes.  Although 
the  method  worked  reasonably  well  on  the  floes  tested,  it  did  not  seem 
as  likely  to  succeed  as  the  feature  tracking.  Hence,  it  was  not  pursued. 

D.  Checking  Results  by  Image  Subtraction 

A  very  useful  technique  was  developed  to  check  the  results  of  ice 
tracking  algorithms.  The  idea  is  first  to  estimate  the  velocity  field  from 
a  given  image  pair.  Then,  using  this  velocity  field,  transport  the  pixels 
in  the  first  image  to  a  new  position  and  form  an  estimate  of  the  second 
image  in  the  pair.  This  second  image  estimate  is  then  subtracted  from 
the  observed  second  image  and  a  difference  image  constructed.  One  can 
then  easily  see  where  the  velocity  field  estimate  is  in  error.  Often  the 
regions  where  there  are  groups  of  errors  are  regions  of  significant  ice 
shear  and  distortion  because  the  ice  is  so  broken  up  that  there  are  no 
coherent  pieces  as  large  as  ten  or  twenty  pixels.  An  example  of  this 
type  of  check  is  given  in  Vesecky  et  al.  ( 1988b). 


IV.  Applications  of  Results 

The  results  of  the  work  performed  here  have  been  primarily  used 
in  the  algorithms  in  place  in  the  ice  tracking  portion  of  the  Geophysical 
Processing  System  (GPS)  at  the  Alaskan  SAR  Facility,  Geophysical 
Institute,  University  of  Alaska,  Fairbanks.  In  particular  we  were  the 
first  to  develop  automated  feature  tracking  and  the  first  to  formulate 
the  basics  of  the  processor  strategy  that  uses  both  feature  tracking  and 
image  pyramid  area  correlation  (IPAC)  techniques.  The  data  products 
from  the  GPS  go  to  Navy  users  via  the  Joint  Ice  Center  (JIC).  The 
primary  data  products  from  the  GPS  are  ice  classification  maps  from 
ERS-1  SAR  images  and  ice  motion  maps  from  pairs  of  ERS-1  SAR 
images.  Data  from  the  Japanese  JERS-1  SAR  are  also  now  becoming 
available. 

The  results  of  our  work  have  also  been  communicated  to  research 
workers  in  the  European  community  for  use  in  analysis  and 
interpretation  or  ERS-1  SAR  images  collected  at  Canadian,  British, 
Norwegian,  German  and  Swedish  ground  stations. 
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V.  Conclusions 

The  primary  conclusions  of  this  research  can  be  summarized  as 
follows: 

1.  Automated  processing  schemes  are  capable  of  estimating  sea  ice 
movement  from  pairs  of  SAR  images  under  a  wide  variety  of  conditions. 

2.  The  image  pyramid  area  correlation  algorithm  when  combined 
with  feature  tracking  methods  can  form  a  robust  ice  velocity  algorithm. 

3.  The  tracking  of  rotating  ice  can  be  handled  by  either  feature 
tracking  or  rotation  invariant  correlation. 

4.  An  excellent  scheme  for  checking  ice  velocity  estimates  is  to  use 
them  to  estimate  the  second  image  of  a  pair  by  transporting  the  pixels 
from  the  first  image.  The  estimate  of  the  second  image  is  then 
subtracted  from  the  observed  second  image  and  regions  of  error  are 
strikingly  displayed. 
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Abstract— Synthetic  aperture  radar  (SAR)  provides  an  excellent 
means  of  observing  the  movement  and  distortion  of  sea  ice  over  large 
temporal  and  spatial  scales.  SAR  observations  are  not  affected  by 
clouds  or  lack  of  sunlight.  The  European  Space  Agency’s  ERS-1  sat¬ 
ellite  as  well  as  others  planned  to  follow  will  carry  SAR’s  over  the  polar 
regions  beginning  in  1989.  A  key  component  in  utilization  of  these  SAR 
data  is  an  automated  scheme  for  extracting  the  sea-ice  velocity  field 
from  a  time  sequence  of  SAR  images  of  the  same  geographical  region. 
We  briefly  describe  two  techniques  for  automated  sea-ice  tracking:  im¬ 
age  pyramid  area  correlation  (hierarchical  correlation)  and  feature 
tracking.  Each  of  these  techniques  is  applied  to  a  pair  of  Seasat  SAR 
sea-ice  images.  The  results  compare  well  with  each  other  and  with 
manually  tracked  estimates  of  the  ice  velocity  field.  We  note  the  ad¬ 
vantages  and  disadvantages  of  each  of  these  automated  methods.  Using 
these  ice  velocity  field  estimates  one  can  reconstruct  one  sea-ice  image 
from  the  other  member  of  the  pair.  Comparing  the  reconstructed  im¬ 
age  with  the  observed  one,  errors  in  the  estimated  velocity  field  can  be 
recognized  and  a  useful  probable  error  display  created  automatically 
to  accompany  ice  velocity  estimates.  We  also  suggest  that  this  error 
display  may  be  useful  in  segmenting  the  sea  ice  observed  into  regions 
that  move  as  rigid  plates  and  regions  of  significant  ice  velocity  shear 
and  distortion. 


I.  Introduction 

HE  IMPORTANCE  of  sea-ice  remote  sensing  arises 
from  the  high  spatial  and  temporal  variability  of  sea 
ice,  its  wide  geographic  extent,  and  its  impact  on  climate, 
arctic  offshore  engineering,  transportation,  and  military 
operations.  Sea  ice  dramatically  changes  the  air-sea  en¬ 
ergy  balance  and  energy  transport.  Since  sea  ice  can  cover 
significant  areas  of  the  Arctic  and  Southern  Oceans,  it  has 
an  important  impact  on  both  polar  and  global  climate. 

The  high  variability  of  sea  ice  coupled  with  its  large 
extent  makes  remote-sensing  techniques  a  necessity  for 
large-scale  studies.  Since  polar  regions  are  often  cloud 
covered  and  without  sunlight  for  long  periods,  the  ability 
of  synthetic  aperture  radar  (SAR)  to  provide  areally  ex¬ 
tensive  high-resolution  images  through  clouds  and  with¬ 
out  sunlight  gives  it  a  central  role  in  sea-ice  remote  sens¬ 
ing.  SAR  provides  good  information  on  sea-ice  extent, 
movement  and  deformation,  internal  geometry  (structure 
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of  floes  and  leads),  and  surface  roughness  [1].  Some  use¬ 
ful  information  on  sea-ice  types  and  their  physical  prop¬ 
erties  can  also  be  obtained  from  SAR  images. 

The  process  of  retrieving  motion  information  from  a 
sequence  of  SAR  images  (of  the  same  geographic  area) 
revolves  around  two  principal  items,  namely: 

1)  Identifying  regions  in  a  sea-ice  field  to  use  as  “tie 
points”  for  ice  floe  tracking. 

2)  Finding  the  same  tie  points  in  subsequent  images  of 
the  sequence. 

Once  the  ice  displacement  is  known  the  ice  velocity  can 
be  calculated  easily  since  the  elapsed  time  between  the 
SAR  images  is  known.  In  this  paper  we  shall  consider  the 
ice  velocity  and  displacement  fields  as  equivalent  since 
they  are  simply  related  by  a  scale  factor. 

This  process  has  been  carried  out  manually  on  Seasat 
SAR  data  by  Thorndike  and  Rothrock  [1]  using  a  split¬ 
screen  video  display  and  an  experienced  operator.  Ice  ve¬ 
locity  is  computed  automatically  once  the  locations  of  tie 
points  are  marked  in  a  pair  of  images  on  the  screen.  Fur¬ 
ther  work  on  satellite  imagery  using  manual  techniques 
has  been  done  by  Carsey  et  al.  [2]  and  by  Curlander,  Holt, 
and  Hussey  [3].  However,  this  process  is  too  slow  and 
expensive  for  real  time  and/or  extensive  analysis  of  im¬ 
ages  provided  by  long-term  satellite  SAR’s.  Two  and  pos¬ 
sibly  three  such  satellite  SAR’s  are  scheduled  to  be  in  or¬ 
bit  by  the  1990’s,  namely  the  European  Space  Agency 
ERS-1,  the  Japanese  ERS-1,  and  the  joint  Canadian-US- 
UK  Radarsat.  Real  aperture  radar  images  at  low  resolu¬ 
tion  are  apparently  available  from  the  Soviet  Kosmos  1500 
series  satellites.  Hence,  the  development  of  automated 
techniques  for  floe  tracking  is  a  primary  requirement  for 
successful  use  of  satellite  SAR  images  in  sea-ice  remote 
sensing. 

This  paper  compares  two  approaches  that  offer  good 
prospects  for  automatically  retrieving  sea-ice  velocity  in¬ 
formation  from  SAR  images  and  presents  a  method  for 
assessing  the  errors  in  the  derived  ice  velocity  field.  The 
two  techniques  compared  here  are  called  “image  pyramid 
area  correlation”  (IPAC),  discussed  in  Section  II,  and 
“boundary  segment  feature  tracking”  (BSFT),  discussed 
in  Section  III.  The  first  technique,  also  known  as  “hier¬ 
archical  correlation,”  comes  from  the  field  of  pattern  rec¬ 
ognition  and  was  first  applied  to  the  tracking  of  sea  ice  by 
Fily  and  Rothrock  [4].  Cross-correlation  of  small  blocks 
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of  pixels  from  a  pair  of  cospatial  images  is  used  to  find 
matching  scene  points.  This  is  made  computationally  ef¬ 
ficient  by  first  building  from  each  primary  high-resolution 
image  a  sequence  (pyramid)  of  derived  images  of  pro¬ 
gressively  coarser  resolution,  but  correspondingly  fewer 
data  points.  Block  cross  correlation  in  the  coarsest  reso¬ 
lution  images  of  the  pyramid  is  computationally  fast  and 
builds  a  crude  model  of  the  displacement  field  of  the  sea 
ice.  This  crude  model  is  then  used  to  guide  searches  in 
the  pyramid  of  progressively  higher  resolution  images.  As 
one  moves  down  the  pyramid  to  the  highest  resolution  im¬ 
age,  searches  are  guided  from  previous  results  and  the 
displacement  field  is  continually  refined.  In  this  paper  we 
discuss  methods  of  rejecting  low-confidence  ice  displace¬ 
ment  estimates  and  interpolating  to  fill  in  the  rejected  data. 
We  also  present  ice  motion  characteristics,  divergence, 
and  curl,  derived  from  the  velocity  field. 

The  second  technique,  discussed  in  Section  III,  also 
comes  from  the  field  of  pattern  recognition.  Our  applica¬ 
tion  to  the  ice  tracing  problem  is  new.  This  method  cor¬ 
responds  more  closely  to  the  manual  (visual)  technique 
than  does  the  image  pyramid  scheme.  Distinctive  features 
are  located  in  the  reference  (initial)  image  and  then  lo¬ 
cated  again  in  the  test  (subsequent)  image.  For  our  fea¬ 
tures  we  use  segments  of  floe-lead  boundaries.  Thus,  each 
image  becomes  a  set  of  boundary  segment  features.  The 
set  of  features  in  the  reference  image  forms  the  tie  points, 
which  are  then  located  in  the  test  image  and  used  to  derive 
the  ice  displacement  (velocity)  field.  We  report  the  results 
of  applying  this  technique  to  Seasat  SAR  images  of  sea 
ice. 

In  Section  IV  we  compare  the  techniques  of  Sections  II 
and  III,  citing  the  advantages  and  disadvantages  of  each. 
Section  V  contains  new  results  applying  image  recon¬ 
struction  and  image  subtraction  to  SAR  sea-ice  image 
pairs.  These  methods  are  used  to  identify  regions  where 
ice  displacement  estimates  are  in  error.  We  also  note  that 
the  regions  of  high  error  density  correspond  to  regions  of 
shear  and  deformation  in  the  sea  ice.  This  paper  ends  with 
Section  VI,  a  summary  of  conclusions  drawn  from  the  re¬ 
sults  reported. 

II.  Image  Pyramid  Area  Correlation  (IPAC) 
Method 

A.  Description  of  Method 

The  problem  of  tracing  sea  ice  from  two  radar  images 
can  be  solved  by  finding  a  correspondence  between  the 
two  images  at  points  that  represent  the  same  scene  point. 
Examples  of  SAR  images  of  sea  ice  are  given  in  Fig.  1, 
where  the  upper  image  is  from  Seasat  orbit  1439  (Oct.  5, 
1978)  and  the  lower  is  from  orbit  1482  (Oct.  8,  1978), 
and  by  Weller  et  al.  [1].  Cross-correlation  of  blocks  of 
pixels  in  the  images  can  be  used  to  find  matching  scene 
points  in  the  images.  The  displacement  vector  between 
matching  points  is  taken  to  be  between  centroids  of  the 
blocks  with  the  largest  correlation  coefficient.  This  dis¬ 
placement  information,  together  with  the  times  of  acqui¬ 
sition  of  the  images,  can  be  used  to  derive  velocity  infor¬ 


mation.  However,  correlation  is  one  of  the  most 
computationally  expensive  image-processing  algorithms. 
One  way  to  speed  up  the  computation  is  to  apply  the  cor¬ 
relation  to  an  image  pyramid  data  structure.  This  was  first 
done  for  sea  ice  by  Fily  and  Rothrock  [4]. 

An  image  pyramid,  shown  in  Fig.  1,  can  be  viewed  as 
a  set  of  images  in  a  stack.  At  the  bottom  of  the  stack  is 
the  original  image  and  at  each  of  the  following  levels  of 
the  stack  is  a  lower  resolution  image  derived  from  the  level 
below  it  by  first  applying  a  low  pass  filter  and  then  resam¬ 
pling  the  image.  The  result  is  a  pyramid  in  which  as  one 
traverses  from  bottom  to  top,  the  number  of  pixels  in  the 
images  decreases  geometrically. 

The  following  steps  are  needed  for  effective  correlation 
in  a  pyramid  structure.  To  begin,  pyramids  (1  and  2)  are 
constructed  (see  Fig.  1);  one  for  each  of  the  two  images. 
Then  the  following  algorithm  is  used: 

1)  For  each  I  x  I  pixel  sample  block  in  the  coarsest 
level  of  pyramid  1,  calculate  a  displacement  vector  by 
finding  the  global  maximum  of  the  correlation  with  the 
coarsest  level  of  pyramid  2.  Generate  a  confidence  mea¬ 
sure  for  the  displacements. 

2)  For  each  finer  resolution  level  of  the  pyramid  do  the 
following: 

a)  Smooth  the  displacement  vector  field  for  the  pre¬ 
vious  (coarser)  pyramid  level  in  order  to  remove 
gross  errors.  This  helps  prevent  errors  in  the  top 
levels  of  the  pyramid  propagating  to  the  lower 
levels. 

b)  Interpolate  the  displacement  vector  field  calcu¬ 
lated  at  the  previous  pyramid  level  to  provide  dis¬ 
placement  estimates  for  each  I  x  I  pixel  sample 
block  of  the  current  level  of  pyramid  1. 

c)  For  each  I  x  I  pixel  sample  block  of  the  current 
level  of  pyramid  1  find  the  local  maximum  of  the 
correlation  with  the  current  level  of  pyramid  2  by 
searching  in  an  M  x  M  window  centered  at  the 
displacement  estimate  of  the  current  level  of  pyr¬ 
amid  2. 

The  algorithm  starts  at  the  top,  coarsest  pyramid  level 
and  proceeds  “top-down”  one  level  at  a  time.  At  the 
coarsest  level,  the  area  to  be  searched  in  the  second  image 
is  not  constrained,  i.e.,  M  is  the  dimension  of  the  entire 
image.  At  each  successive  level  in  the  pyramid,  the  esti¬ 
mated  location  of  the  largest  correlation  found  at  the  pre¬ 
vious,  coarser  pyramid  level  constrains  the  locations 
searched.  These  location  constraints  greatly  decrease  the 
total  number  of  locations  that  must  be  searched.  The  over¬ 
all  gain  in  speed  over  brute-force  correlation  is  between 
100  and  1000. 

The  selection  of  the  size  (I)  of  the  correlation  sample 
block  involves  a  trade-off  between  accurate  determination 
of  the  location  of  correlation  maximum  (small  I)  and 
avoidance  of  false  correlation  maxima  (large  I).  Vesecky 
et  al.  [5]  describe  a  cross-correlation  technique  for  choos¬ 
ing  I.  The  choice  of  I  is  dependent  on  image  content.  For 
the  sea-ice  image  of  Fig.  1 ,  an  I  corresponding  to  about 
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Fig.  !.  Image  pyramids  for  two  Seasat  SAR  images  of  the  same  central 

pack  ice  area  in  the  Beaufort  Sea.  The  images  were  collected  three  days 

apart  on  Oct.  5  and  8,  1978  on  orbits  1439  and  1482. 

a  4  x  4  km  block  size  was  judged  to  be  near  optimum, 
i.e.,  I  =  16  for  250-m  size  pixels.  This  size  is  of  the  order 
of  the  size  of  the  average  ice  floe. 

The  quantity  M  is  selected  to  correspond  to  a  search 
area  of  the  order  of  the  physical  size  of  the  correlation 
sample  block  in  the  coarser  pyramid  level  just  processed. 
We  have  taken  M  such  that  the  search  area  is  four  times 
larger  then  the  previous  (coarser)  level  sample  block  size. 

B .  First-Order  Results  for  Seasat  SAR  Image  Pair 

We  applied  the  pyramid  algorithm  described  above  to 
two  averaged  (250-m  pixels)  Seasat  images  of  central  arc¬ 
tic  pack  sea  ice-Rev.  1439  (October  5,  1978)  and  Rev. 
1482  (October  8,  1978)  [6].  The  experiment  used  three- 
level  pyramids.  A  median  filter  was  used  for  the  smooth¬ 
ing  of  the  displacement  vectors.  Bilinear  interpolation  was 
used  for  the  interpolation  of  the  displacement  vectors.  Fig. 
2(a)  shows  the  displacement  vectors  at  the  finest  resolu¬ 
tion.  The  direction  of  most  vectors  is  from  top  to  bottom. 
Comparing  these  displacement  vectors  with  a  computer- 
assisted  manual  analysis  by  Thorndike  and  Rothrock  [1] 
shown  in  Fig.  2(b),  it  is  seen  that  the  method  performs 
very  well.  However,  there  are  areas  where  the  automated 
algorithm  has  generated  erroneous  velocity  vectors.  The 
problem  areas  at  the  edges  of  the  figure  correspond  to  ice 
features  flowing  into  or  out  of  the  finite  area  of  these  im¬ 
ages. 

C.  Confidence  Estimates 

There  are  two  parts  of  the  pyramid  algorithm  where 
confidence  measures  are  useful.  For  the  final  vector  out¬ 
put,  confidence  measures  can  be  used  to  prune  or  to  at 
least  highlight  suspect  vectors.  In  addition,  during  the  in¬ 
terpolations  from  a  coarser  pyramid  level  to  a  finer  pyra¬ 
mid  level,  the  confidence  measures  can  be  used  as  a 


(b) 

Fig.  2.  (a)  Sea-ice  velocity  vectors  generated  by  the  image  pyramid  area 
correlation  technique  (IPAC)  applied  to  the  SAR  images  of  Fig.  1.  (b) 
Sea-ice  velocity  vectors  generated  from  the  same  image  data  of  Fig.  1 
by  Thorndike  and  Rothrock  {1]  using  a  computer-assisted  manual  tech¬ 
nique.  The  velocity  estimates  of  (a)  are  improved  by  pruning  of  low 
confidence  estimates  and  replacement  by  interpolation  as  shown  in  Fig. 
4. 

weighting  factor  for  the  interpolation  [7].  The  first  tech¬ 
nique  was  used  in  obtaining  the  results  reported  here. 

The  coordinates  of  the  maximum  correlation  for  a  block 
provide  the  value  for  its  velocity  vector.  The  values  of  the 
correlation  surface  near  the  maximum  may  provide  a  con¬ 
fidence  measure  to  be  associated  with  the  velocity  vec¬ 
tors.  Fig.  3  shows  the  value  of  the  correlation  coefficient 
superimposed  on  the  velocity  vectors.  The  correlation 
coefficient  is  being  used  as  a  measure  of  confidence,  with 
brighter  regions  signifying  high  confidence  and  darker  re¬ 
gions  signifying  low  confidence.  It  is  seen  that  there  is  a 
correspondence  between  areas  of  low  confidence  and  er¬ 
rors  in  the  displacement  field. 
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Fig.  3.  Confidence  estimates  for  sea-ice  velocity  vectors  generated  by  the 
automated  image  pyramid  area  correlation  technique  (IPAC).  The  bright 
shades  denote  regions  of  high  confidence  and  the  dark  areas,  regions  of 
low  confidence.  The  ice  velocity  field  is  superimposed  on  the  confidence 
estimates,  which  are  based  on  the  magnitude  of  the  cross-correlation  coef¬ 
ficient  for  that  particular  velocity  estimate. 


D.  Smoothing  Techniques 

In  the  “real"  world  of  sea-ice  remote  sensing  using 
SAR  images  one  encounters  noise,  grey  level  variation 
across  a  given  image  and  between  images,  and  gross  geo¬ 
metric  distortions.  These  problems  result  in  errors  when 
block  correlation  is  used  for  finding  ice  displacement  vec¬ 
tors.  In  view  of  these  problems,  the  assumption  is  usually 
made  that  the  displacement  field  is  smooth  for  most  areas 
of  the  image  [8].  Smoothing  techniques  may  be  used  to 
remove  the  gross  errors  in  the  displacement  field.  These 
techniques  are  especially  important  for  pyramid  correla¬ 
tion  because  errors  made  at  the  coarser  resolutions  prop¬ 
agate  to  the  finer  levels  due  to  the  constraints  on  search 
locations  at  the  finer  levels.  One  must  be  careful  when 
using  smoothing  techniques,  however,  in  areas  where 
there  is  high  shear. 

In  the  algorithm  described  in  Section  II-A  smoothing  of 
the  estimated  ice  velocity  field  was  done  at  each  pyramid 
level  by  the  following  procedure.  The  coarse  level  vectors 
were  replaced  by  the  median  of  the  surrounding  high  con¬ 
fidence  vectors.  Bilinear  interpolation  was  then  used  to  fill 
in  vectors  for  the  next  higher  resolution  level  in  the  pyr¬ 
amid.  At  the  highest  resolution  level  of  the  pyramid  this 
procedure  produced  the  result  of  Fig.  2(a). 

To  improve  the  result  of  Fig.  2(a)  we  have  implemented 
a  different  smoothing  technique  by  extending  the  con- 
trolled-continuity  splines  method  of  Terzopoulos  [9]  for 
fitting  vector  fields.  This  more  versatile  interpolation 
technique  makes  systematic  use  of  the  confidence  mea¬ 
sure  and  minimizes  the  weighted  square  deviation  from 
the  displacement  data  while  at  the  same  time  minimizing 
the  integral  of  the  second  spatial  derivatives.  The  solution 


Fig.  4.  Icc  velocity  field  of  Fig.  2(a)  after  pruning  the  low  confidence  vec¬ 
tors  and  using  the  control led-continuity  spline  technique  for  interpola¬ 
tion. 


to  the  minimization  problem  is  via  finite  elements.  The 
resulting  matrix  equation  is  solved  via  a  multigrid  relax¬ 
ation  technique  using  the  Gauss-Seidel  method  for  intra¬ 
level  computations. 

Fig.  4  shows  the  velocity  field  after  pruning  the  low 
confidence  vectors  and  using  controlled-continuity  splines 
for  interpolation.  This  surface  fitting  technique  could  be 
extended  to  allow  discontinuities  in  the  smoothing  pro¬ 
cess  and  thus  provide  a  viable  smoothing  method  for  areas 
with  velocity  shear. 

E .  Divergence  and  Curl 

We  have  evaluated  the  divergence  and  curl  of  the  com¬ 
puted  displacement  field,  which  may  be  considered  to  be 
an  averaged  velocity  field.  The  displacement  field  was  first 
smoothed  as  in  Fig.  4;  numerical  derivatives  were  then 
computed.  This  smoothing  reduces  the  effect  of  finite 
block  size  on  estimates  of  spatial  derivatives.  Fig.  5(a) 
shows  the  divergence  of  the  velocity  field.  Sinks  are 
shown  as  dark  areas  in  the  image  and  sources  are  shown 
as  bright  areas  in  the  image.  Within  rigid  blocks,  the  in¬ 
termediate  gray  level  signifies  zero  divergence,  as  ex¬ 
pected.  At  block  boundaries,  however,  the  relative  mo¬ 
tions  of  the  blocks  result  in  nonzero  divergence. 
Comparing  the  divergence  plot  of  5(a)  with  the  sea-ice 
images  of  Fig.  1  we  note  that  the  bright  source  region  in 
the  lower  left  quadrant  corresponds  to  a  region  where  the 
open  water  area  increases  between  the  5  (upper)  and  8 
(lower)  October  images.  A  similar  source  feature  occurs 
at  the  upper  center  of  5(a).  Fig.  5(b)  shows  the  curl  of  the 
field.  Since  the  curl  of  a  two-dimensional  field  has  only 
one  component,  it  can  be  represented  as  a  scalar.  The  curl 
at  the  center  of  a  clockwise  rotation  is  shown  as  a  brighter 
area.  The  curl  at  the  center  of  a  counterclockwise  rotation 
is  shown  as  a  darker  area. 


42 


Ifchfc  TRANSACTIONS  ON  GEOSCIENCE  AND  REMOTE  SENSING.  VOL  26.  NO  I.  JANUARY  1988 


<b) 


Fig.  5.  Divergence  (a)  and  curl  (b>  of  the  smoothed  ice  velocity  field  of 
Fig.  4.  In  (a)  sinks  are  shown  as  dark  areas  in  the  image  and  sources  as 
bright  areas.  Gray  areas  signify  zero  divergence.  Since  the  curl  of  a  two- 
dimensional  field  has  only  one  component,  it  can  be  represented  as  a 
scalar.  The  curl  at  the  center  of  a  clockwise  rotation  is  shown  as  a  brighter 
area.  The  curl  at  the  center  of  a  counterclockwise  rotation  is  shown  as  a 
darker  area. 

III.  Feature  Tracking  With  Floe-Lead  Boundaries 
A.  Introduction  to  the  Method 
The  second  algorithm  approaches  the  problem  of  sea- 
ice  tracking  somewhat  differently  by  relying  on  selected 
sea-ice  features  rather  than  on  the  content  of  the  entire 
image.  There  are  a  wide  variety  of  features  in  sea  ice  that 
one  could  track,  including  pressure  ridges,  isolated  floes, 
and  floe-lead  boundaries.  Here  we  have  chosen  floe-lead 
boundaries  as  the  physical  feature  we  will  track  because 
they  are  usually  the  most  prominent,  distinct,  and  high 
constrast  features  in  SAR  images  of  sea  ice.  The  move¬ 
ment  of  these  features  from  one  image  to  the  next  pro¬ 
vides  the  information  from  which  the  sea-ice  velocity  field 


is  estimated.  Feature  patterns  must  be  very  distinctive,  if 
not  unique,  objects.  They  are  characterized  by  informa¬ 
tion  that  is  independent  of  the  object  location  and  orien¬ 
tation.  though  this  latter  information  is,  of  course,  at¬ 
tached  to  them. 

In  the  feature  tracking  method  each  image  is  trans¬ 
formed  into  a  set  of  features  described  by  parameters  in¬ 
dependent  of  the  feature  location  and  orientation.  The  col¬ 
lection  of  several  hundred  features,  each  described  by 
several  hundred  numbers,  reduces  the  information  needed 
to  characterize  a  100  x  100  km  image  by  roughly  two 
orders  of  magnitude.  The  search  space  is  correspondingly 
reduced. 

Once  a  feature  pattern  is  recognized  in  both  members 
of  a  pair  of  images,  sea-ice  displacement  and,  hence,  ve¬ 
locity  information,  can  be  derived.  Since  we  know  the 
movement  of  every  pixel  in  the  feature,  we  can  derive  not 
just  one,  but  many  velocity  vectors.  Further,  we  know  the 
rotation  as  well  as  the  displacement  of  a  feature,  since  we 
keep  track  of  a  feature’s  geographic  orientation  in  both 
images  of  the  pair.  The  results  shown  here  are  for  floe- 
lead  boundaries;  however,  the  method  can  also  be  applied 
to  other  features  such  as  pressure  ridges  or  isolated  floes. 

B.  Description  of  the  Feature  Tracking  Algorithm 

The  feature  tracking  algorithm  can  best  be  understood 
by  first  considering  the  overview  in  Fig.  6  and  then  dis¬ 
cussing  each  component  in  turn.  Following  Fig.  6.  each 
image  is  first  averaged  to  a  resolution  of  about  250  m.  Ice 
and  lead  are  then  classified  by  a  simple  intensity  threshold 
and  floe-lead  boundaries  are  identified.  Pixels  along  floe- 
lead  boundaries  are  classified  and  assembled  into  seg¬ 
ments.  These  boundary  segments  are  characterized  and 
thus  form  the  features  that  are  tracked  from  one  SAR  im¬ 
age  to  another  collected  three  days  later.  Several  tracking 
passes  are  made  through  the  collections  of  features  in  each 
image.  The  high  confidence  ice  displacement  estimates, 
corresponding  to  high  correlation  between  the  character¬ 
istics  of  matching  features  in  the  two  images,  are  used  to 
guide  the  search  for  other  matching  features  and  thus  more 
ice  displacement  estimates. 

Referring  to  the  algorithm  block  diagram  of  Fig.  6,  we 
first  average  both  SAR  images  of  a  pair  to  increase  the 
pixel  size  from  the  original  25  x  25  m  size  to  250  x  250 
m  size.  Averaging  reduces  the  speckle  variance  inherent 
in  SAR  images  and  make  computation  more  tractable  since 
there  are  a  factor  of  100  fewer  pixels  to  deal  with  after 
averaging.  Since  SAR  pixels  cannot  be  located  geograph¬ 
ically  to  better  than  about  200  m,  the  average  process  re¬ 
sults  in  little  loss  of  velocity  measurement  accuracy. 

Proceeding  down  the  block  diagram  of  Fig.  6,  a  histo¬ 
gram  of  pixel  intensity  is  formed  and  a  simple  pixel  in¬ 
tensity  threshold  is  chosen  to  distinguish  between  ice  floes 
and  leads.  Fily  and  Rothrock  [4]  have  shown  that  this 
technique  is  viable  for  an  example  SAR  image  of  central 
pack  ice.  Thus,  the  SAR  image  pair  can  be  transformed 
into  a  pair  of  images  with  binary  pixels,  i.e.,  1  for  ice  floe 
pixels  and  0  for  lead  water  pixels  (including  new  ice).  In 
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Fig.  6.  Block  diagram  of  feature  tracking  method  using  hoe-lead  bound¬ 
aries. 

the  next  block  of  Fig.  6  we  identify  floe-lead  boundary 
pixels  simply  by  finding  water  pixels  adjacent  to  ice  pix¬ 
els. 

In  the  next  block  of  Fig.  6  the  ice-water  boundary  pix¬ 
els  identified  above  are  classified  into  eight  types  accord¬ 
ing  to  their  surroundings  to  determine  the  orientation  of 
the  floe-lead  boundary  edge.  The  four  nearest  neighbor 
pixels  are  considered  to  find  if  there  is  water  above,  be¬ 
low,  or  to  the  right  or  to  the  left  of  a  given  boundary  pixel. 
This  scheme  is  illustrated  in  Fig.  7. 

The  classified  pixels  are  then  assembled  into  boundary 
segments  that  are  the  features  to  be  tracked.  The  assembly 
scheme  used  here  begins  by  finding  the  dominant  bound¬ 
ary  direction.  This  is  done  by  finding  the  dominant  edge 
pixel  types  of  Fig.  7.  For  example,  if  types  3  and  7  are 
most  common  the  dominant  boundary  direction  is  vertical 
in  Fig.  7.  In  this  example  we  would  then  divide  our 
boundary  segments  into  two  groups,  corresponding  to  the 
left  and  right  edges  of  the  leads  in  Fig.  1,  to  wit: 

Left  Edges:  edge  pixel  types  2,  3,  4  and  associated 
types  1  and  5 

Right  edges:  edge  pixel  types  6,  7,  8  and  associated 
types  1  and  5. 

Our  aim  then  would  be  to  match  left  edges  in  one  image 
with  left  edges  in  the  other  image  and.  independently, 
right  edges  in  the  two  images.  This  allows  the  ice  floes 
along  the  left  side  of  a  lead  to  slip  with  respect  to  ice  floes 
along  the  right  side.  Thus,  shear  motion,  which  we  would 


Fig.  7.  Diagram  illustrating  edge  pixel  classification  scheme.  The  black 
circles  represent  water  pixels  and  the  white  circles  ice  pixels. 

expect  across  leads  in  central  pack  ice,  is  built  into  the 
algorithm  for  this  type  of  sea  ice. 

In  the  next  block  of  Fig.  6  boundary  segment  features 
are  characterized  by  their  deviation  from  a  straight  line. 
First,  a  "‘least  squares"  criterion  straight  line  is  fit  to  a 
given  boundary  segment  assembled  as  above.  The  dis¬ 
tance  from  this  straight  line  to  the  segment  is  calculated 
at  regular  intervals  along  the  straight  line.  This  set  of  de¬ 
viations  characterizes  the  shape  of  a  boundary  segment  as 
illustrated  in  Fig.  8.  Note  that  this  method  of  character¬ 
izing  a  boundary  segment  feature  is  invariant  under  rota¬ 
tion  of  the  feature. 

Recognition  of  matching  features  (Fig.  6)  is  done  sim¬ 
ply  by  comparing  the  sets  of  deviation  measurements  for 
two  candidate  features.  A  variance  normalized  cross-cor¬ 
relation  coefficient  ( C )  was  used  to  evaluate  how  closely 
two  candidate  features  matched.  Often  segments  that  we 
want  to  compare  are  of  different  lengths  because  of  partial 
breakup  of  a  boundary  segment  over  the  three  days  be¬ 
tween  images.  So  in  our  search  we  include  features  of 
different  lengths.  During  the  comparison  process  we  slide 
the  shorter  feature  along  the  longer  feature  with  C  calcu¬ 
lated  for  all  placements  of  the  smaller  feature  with  respect 
to  the  larger.  This  allows  for  larger  segments  to  break  into 
smaller  ones  and  still  be  found. 

We  have  at  this  point  (the  entrance  to  the  loop  structure 
in  Fig.  6)  two  images  that  each  contain  sets  of  boundary 
segments.  Because  in  the  central  pack  ice  we  do  not  ex¬ 
pect  large  rotations  we  have  divided  each  of  these  sets 
into  left  and  right  edges.  The  search  procedure  must  select 
a  feature  from  the  reference  image  (first  in  time)  and  try 
to  match  it  to  a  feature  in  the  test  image  (later  in  time).  A 
successful  match  reveals  ice  displacement  and,  hence,  ice 
velocities  at  a  number  of  points  all  along  the  matching 
segments. 

The  search  and  matching  procedure  (loop  structure  in 
Fig.  6)  is  an  iterative  one  and  has  features  that  remind  one 
of  the  image  pyramid  method  of  Section  II  above.  During 
the  several  iterations  of  the  search  we  use  the  high  con¬ 
fidence  ice  displacement  estimates  of  the  initial  iterations 
to  guide  the  searches  of  subsequent  iterations.  In  the  first 
iteration  only  long  boundary  segments  are  used  since  the 
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Fig.  8.  Boundary  segment  fitted  by  a  straight  lire  and  characterized  b.* 
deviation  from  the  fitted  line. 


larger  number  of  pixels  gives  a  more  stable  estimate  of  C 
and  a  correct  match  yields  ice  displacement  estimates  over 
a  larger  region.  At  this  stage  for  a  left  boundary  segment 
of  n  pixels  in  the  reference  image  we  search  over  left 
boundary  segments  with  n/2  to  2 n  pixels.  The  search  is 
limited  to  a  displacement  of  60-20  km  per  day  for  the 
three  days  between  the  images.  This  search  area  could  be 
significantly  reduced  if  environmental  data  such  as  local 
winds  were  available  to  indicate  probable  ice  movement. 

In  the  first  iteration  we  want  to  have  a  very  high  con¬ 
fidence  in  our  ice  displacement  estimates.  This  is  because 
we  will  use  them  to  guide  searches  in  subsequent  itera¬ 
tions.  Hence,  we  require  high  values  of  C  for  a  match  and 
reject  matches  that  are  not  confirmed  by  nearby  matches 
having  roughly  similar  displacements,  +  /  -20°  in  direc¬ 
tion  and  +  /  -  30  percent  in  magnitude.  In  subsequent  it¬ 
erations  smaller  segments  are  used,  the  search  region  is 
smaller  and  guided  by  previous  results,  and  lower  values 
of  C  are  demanded  for  a  match.  In  the  results  reported 
here  seven  iterations  were  made. 

C.  Application  to  Seasat  SAR  Images  in  the  Beaufort 
Sea 

We  now  apply  the  boundary  segment  feature  tracking 
algorithm  to  the  same  set  of  sea-ice  images  in  Fig.  1 .  We 
use  only  the  250  x  250  m  pixel  size  images.  Boundary 
segments  constructed  from  these  images  vary  in  length 
from  4  to  163  pixels  or  1  to  41  km  in  length  and  have  a 
dominant  direction  which  is  vertical,  i.e.,  North-South, 
in  Fig.  1.  The  resulting  sets  of  boundary  segment  features 
in  the  reference  and  test  images  are  shown  in  Fig.  9. 

In  the  initial  iteration  only  the  longer  segments  of  im¬ 
age  1439  with  lengths  greater  than  30  pixels  (7.5  km)  were 
considered.  Thus,  in  the  test  image  (orbit  1482)  matching 
segments  of  length  n/2  =  15  pixels  and  above  were  con¬ 
sidered  in  the  search.  The  rather  strict  criteria  on  match¬ 
ing  features  were  successful  in  removing  many  erroneous 
matches,  although  it  did  also  eliminate  a  few  correct 
matches  as  well.  The  subsequent  iterations  added  pro¬ 
gressively  fewer  matches  until  by  the  seventh  iteration  no 
new  matches  were  found.  It  has  been  our  experience  in 
developing  this  algorithm  that  the  primary  difficulty  lies 
in  making  sure  incorrect  matches  are  rejected  rather  than 
in  making  sure  correct  matches  are  retained. 

Fig.  10  shows  the  resulting  ice  displacement  vectors  for 
the  sea-ice  image  pair  of  Fig.  1.  These  displacement  vec¬ 
tors  agree  very  well  with  the  manually  generated  and 
IPAC  displacement  vectors  of  Fig.  2.  We  note  that  there 
are  large  regions  where  no  displacement  vectors  are 
shown.  This  situation  arises  from  two  causes.  First,  no 
ice  displacement  vectors  will  be  generated  in  regions 
where  there  are  no  lead-floe  boundary  features  to  be 
tracked.  The  second  is  that  in  regions  where  there  is  too 


much  deformation,  a  boundary  segment  in  the  reference 
image  will  often  be  too  badly  broken  up  in  the  test  image 
to  be  recognized.  In  the  first  case  other  features,  such  as 
pressure  ridge  patterns,  could  be  tracked  as  well  as  the 
lead-floe  boundary  segments.  In  the  second  case  the  use 
of  higher  resolution  images  might  be  helpful  if  the  defor¬ 
mation  were  less  severe  at  smaller  size  scales.  A  draw¬ 
back  of  using  higher  resolution  pixels  in  SAR  images  is 
that  speckle  noise  becomes  a  significant  component  of  the 
intensity  variation  from  one  pixel  to  the  next. 

IV.  Comparison  of  Methods 
Given  our  results,  we  can  qualitatively  address  some  of 
the  advantages  and  disadvantages  that  accompany  the  sea- 
ice  tracking  methods  considered  above.  In  particular,  we 
can  compare  the  methods  using  an  outline  suggested  by 

Kashef  and  Sawchuck  [10]. 

A  feature- matching  method,  such  as  the  boundary  seg¬ 
ment  feature  tracking  (BSFT)  method  of  Section  III,  ap¬ 
proaches  mat  (image)  data  differently  from  a  correlation 
method,  such  as  the  image  pyramid  area  correlation 
(IPAC)  method  of  Section  II.  Whereas  correlation  meth¬ 
ods  rely  on  the  intensity  values  of  map  data,  feature¬ 
matching  methods  rely  on  detected  edges,  boundaries,  and 
vertices  of  map  data.  The  advantage  of  using  BSFT  in 
tracking  sea  ice  lies  in  reducing  an  image  of  millions  of 
pixels  to  one  of,  say,  thousands  of  pixels.  Only  relevant 
data  applies  to  tracking.  The  potential  cost-savings  for 
computational  speed  exist,  depending  on  how  “relevant 
the  reduced  data  are.  In  particular,  one  should  only  pick 
those  edges,  boundaries,  and  vertices  common  to  both  im¬ 
ages  for  tracking. 

Picking  those  edges,  boundaries,  and  vertices  common 
to  both  images  remains  difficult.  Thus,  the  advantage  for 
BSFT  also  harbors  a  disadvantage.  Boundaries  may  frag¬ 
ment  or  disappear.  Furthermore,  unrelated  boundaries 
may  have  similar  characteristics  and  lead  to  false  matches. 
Finally,  some  areas  of  sea  ice  are  virtually  featureless, 
having  few  discemable  floe-lead  boundaries. 

To  address  the  difficulty  in  identifying  features  com¬ 
mon  to  both  images,  we  have  adopted  a  feature  tracking 
algorithm  that  iterates.  For  the  first  pass,  feature  tracking 
results  in  finding  five  to  ten  matching  features.  We  use 
these  matches  to  guide  the  search  for  other  matches  on  the 
next  pass.  The  process  repeats.  In  this  way,  we  have  suc¬ 
cessfully  identified  a  significant  number  of  features  com¬ 
mon  to  images  for  orbit  1439  and  orbit  1482.  Iterative 
feature  tracking  has  resulted  in  a  velocity  vector  map  that 
agrees  well  with  manual  tracking  results,  although  the 
number  of  velocity  vectors  located  by  feature  tracking  falls 
short  of  the  number  of  velocity  vectors  located  either  by 
IPAC  or  manual  tracking.  In  die  case  at  hand  we  expect 
this  behavior  since  only  one  type  of  feature  was  consid- 
ered. 

Correlation  methods  surmount  the  difficulty  of  identi- 
fying  common  features  or  tie  points  by  allowing  a  large 
range  of  image  content:  i.e.,  the  algorithms  are  not  de¬ 
pendent  on  the  existence  of  specific  features  such  as  floe- 
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Fig.  9.  Search  space  of  long  floe-lead  boundary  segments  constructed  from  orbit  1439  (left)  and  orbit  1482  tright)  SAR  image  data.  Only  left  edges  are 

shown. 


Fig.  10.  Ice  velocity  vectors  derived  from  the  SAR  images  of  Fig.  I  using 
the  boundary  segment  feature  tracking  (BSFT)  algorithm. 


lead  boundaries.  Traditionally,  computation-intensive 
correlation  methods  become  attractive  when  computation 
is  reduced  by  the  IPAC  approach.  The  noise  immunity 
provided  by  averaging  in  IPAC  comes  as  an  added  bonus. 
Therefore,  the  advantage  of  using  IPAC  techniques  stems 
threefold:  IPAC  assumes  little  about  image  content;  IPAC 
can  track  global  displacements  well;  IPAC  has  inherent 
noise  immunity. 

Despite  the  advantages  of  IPAC,  “image  pyramid" 
techniques,  like  the  one  used  here,  cannot  detect  small- 
scale  pattern  changes  that  involve  a  large  rotation.  This 
problem  is  significant  given  the  rotational  as  well  as  trans¬ 
lational,  nature  of  sea-ice  motion.  To  address  this  draw¬ 
back,  one  can  opt  for  two  schemes.  One  scheme  involves 
substituting  at  some  point  a  rotation-invariant  correlation 


measure  in  place  of  standard  correlation.  Another  ap¬ 
proach  would  be  a  synthesis  from  both  IPAC  and  BSFT. 

There  are  other  disadvantages  to  using  “image  pyra¬ 
mid"  techniques.  In  our  research,  we  have  found  that 
IPAC  was  sensitive  to  tracking  mistakes  at  the  top  (course 
resolution  level)  of  the  pyramid.  These  tracking  errors 
propagated  and  resultecfin  some  velocity  vector  estimates 
significantly  different  than  manual  tracking  results.  For 
the  central  pack  ice  images  we  have  used,  we  have  been 
able  to  discard  a  velocity  vector  if  the  correlation  coeffi¬ 
cient  corresponding  to  it  falls  below  a  threshold  value  and 
then  fill  the  void  by  an  interpolation  scheme.  The  scheme 
of  discarding  velocity  vectors  judged  to  be  unreliable  and 
interpolating  from  surrounding  vectors  to  fill  in  the  void 
can  also  be  applied  to  BSFT  ice  velocity  estimates.  Un¬ 
fortunately,  the  same  method  might  not  work  so  well  for 
other  sea-ice  images;  for  example,  if  voids  occur  over 
areas  too  large  to  interpolate  across.  Other  potential  prob¬ 
lems  using  “image  pyramid"  techniques  might  arise  in 
areas  of  nonuniform  motion  from  averaging  velocity  vec¬ 
tor  estimates  over  too  many  pixels. 

Because  IPAC  and  BSFT  provide  estimates  of  ice  ve¬ 
locity  that  are  independent  the  two  methods  are  comple¬ 
mentary  in  the  sense  that  agreement  between  the  two 
methods  implies  a  high  confidence  estimate  of  ice  veloc¬ 
ity.  Thus,  the  two  methods  could  be  used  together  with  a 
single  algorithm  to  provide  better  sea-ice  velocity  esti¬ 
mates  than  either  could  alone.  For  example,  BSFT  veloc¬ 
ity  estimates  could  be  used  to  reduce  critical  IPAC  errors 
in  the  coarse  resolution  levels  of  an  image  pyramid  and 
to  confirm  results  at  finer  resolution  levels. 

V.  Error  Analysis  by  Image  Reconstruction  and 
Subtraction 

A.  Image  Restoration .  Remapping ,  and  Image 
Subtraction 

As  stated  in  the  introduction,  the  goal  of  this  research 
is  to  produce  ice  velocity  data  for  sea-ice  science  and  ap- 
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plications.  To  make  the  ice  velocity  field  estimates  of 
Figs.  2,  4.  and  10  truly  useful  they  must  be  accompanied 
by  some  measure  of  the  quality  of  the  estimates.  We  have 
found  image  remapping  to  be  a  useful  technique  for  eval¬ 
uating  the  quality  of  a  velocity  field.  The  estimated  ice 
velocity  field  defines  a  mapping  between  the  pair  of  im¬ 
ages.  One  can  use  this  mapping  to  remap  the  intensities 
in  the  second  image  to  the  locations  they  would  be  found 
in  the  first  image.  When  we  alternatively  display  the  first 
image  and  the  remapped  second  image  on  the  same  image 
display,  the  areas  where  the  algorithm  works  well  are  seen 
to  be  static.  The  areas  where  the  algorithm  fails  appear  to 
flicker  or  to  move  slightly  because  of  errors  in  the  map 

ping.  ._ 

With  the  remapped  data  available,  the  absolute  differ¬ 
ence  of  the  first  image  and  the  remapped  second  image 
can  be  produced.  In  order  to  understand  the  significance 
of  the  difference  technique,  the  motion  detection  problem 
must  be  expressed  as  an  image  restoration  problem. 

We  have  modeled  the  motion  process  as  a  linear  space- 
varying  transformation  of  intensities  from  the  first  image 
to  the  second  image.  Using  this  model,  the  second  image 
may  be  considered  to  be  a  superposition  integral  involv¬ 
ing  the  impulse  response  of  the  motion  process  and  the 
intensity  function  of  the  first  image.  We  further  assume 
that  the  imaging  system  does  not  vary  from  image  to  im¬ 
age,  and  that  the  correlation  algorithm  can  also  be  mod¬ 
eled  as  a  space- varying  linear  process.  Then,  it  can  be 
shown  that  the  difference  between  the  first  image  and  the 
remapped  version  of  the  second  image  bounds  the  errors 
in  the  motion  if  there  is  detail  in  the  image.  Intuitively,  if 
the  motion  algorithm  is  correct,  the  remapped  intensities 
subtract  perfectly  from  the  original  intensities  leaving  zero 
remainder.  If  it  is  erroneous,  then  the  intensities  will  not 
subtract  perfectly  and  the  residual  errors  will  be  displayed 
as  bright  portions  in  the  difference  image. 

We  recognize  that  in  the  practical  world  of  SAR  remote 
sensing  there  are  likely  to  be  changes  in  intensity  between 
the  two  images  in  a  pair  that  are  unrelated  to  ice  motion. 
For  example,  changes  in  system  gain,  incidence  angle, 
aspect  angle,  and  surface  conditions  (such  as  snow  and 
wind)  can  cause  variations  in  SAR  image  intensity.  How¬ 
ever,  the  successful  results  of  Fig.  11  resulted  from  the 
pair  of  Seasat  SAR  images  shown  in  Fig.  1.  So  we  are 
optimistic  that  the  method  will  have  application  to  a  sig¬ 
nificant  fraction  of  satellite  SAR  observations  of  sea  ice, 
which  are  less  subject  to  the  variations  mentioned  above 
than  aircraft  observations. 

B.  Error  Detection  by  Image  Subtraction 

Fig.  1 1  shows  the  absolute  value  of  the  difference  be¬ 
tween  the  Oct.  5  (orbit  1439)  image  and  the  remapped 
Oct.  8  (orbit  1482)  image.  The  remapped  Oct.  8  image 
was  created  using  the  smoothed  IP  AC  velocity  vector  es¬ 
timates  of  Fig.  4.  The  difference  picture  has  been  dis¬ 
played  so  that  image  intensity  differences  greater  than  a 
chosen  threshold  value  are  bright  and  those  below  are 


Fig  1 1  Image  reconstruction  and  subtraction  map  based  on  the  ice  veloc¬ 
ity  field  of  Fig.  4.  This  difference  picture  is  displayed  so  that  absolute 
image  intensity  differences  greater  than  a  chosen  threshold  value  are 
bright  and  those  below  are  dark.  Thus,  regions  where  ice  velocity  esti- 
mates  are  likely  to  be  in  error  are  bright  and  low  error  regions  are  dark. 


dark.  Thus,  high  error  regions  are  bright  and  low  error 
regions  are  dark. 

This  image  provides  a  convenient  display  of  the  regions 
where  the  ice  velocity  estimates  are  likely  to  be  in  error. 
Here  we  have  assumed  that  the  ice  itself  remains  the  same, 
as  seen  by  the  SAR.  in  the  two  images.  This  error  map 
can  be  used  in  two  ways.  First,  it  can  be  used  as  a  guide 
to  the  credibility  of  the  ice  velocity  field  estimates  pro¬ 
vided  by  SAR.  Thus,  the  ice  science  or  applications  user 
is  in  a  position  to  know  when  and  where  the  ice  velocity 
data  provided  to  him  are  dependable  or  likely  to  be  in 
error.  Second,  the  error  data  can  be  used  as  feedback 
within  an  ice  velocity  estimation  algorithm  to  improve  the 
velocity  estimates.  For  example,  higher  resolution  image 
data,  i.e,  pixels  smaller  than  the  250  m  used  here,  could 
be  employed  to  try  to  resolve  errors  sensed  by  the  image 
subtraction  technique. 

C.  Interpretation  of  Regions  Containing  Large  and 
Small  Errors 

By  comparing  the  error  display  of  Fig.  1 1  with  the  ice 
velocity  field  as  shown  in  Figs.  2,  4,  and  10  we  find  that 
we  can  make  an  interesting  interpretation  of  the  error  dis¬ 
play.  First,  we  note  that  the  large  dark  areas  where  errors 
are  low  correspond  to  large  blocks  of  ice  that  move  as 
rigid  plates.  These  large  units  were  also  pointed  out  by 
Fily  and  Rothrock  [4].  Second,  we  see  that  the  large  error 
(bright)  regions  tend  to  fall  along  lines  that  separate  the 
large  rigid  plates.  It  is  in  these  regions  that  the  large  plates 
apparently  grind  against  one  another  producing  velocity 
shear  that  deforms  the  ice.  The  deformation  of  the  ice  in 
these  regions,  coupled  with  the  fact  that  the  ice  velocity 
estimates  are  discrete  and  smoothed,  causes  errors  in  ice 
velocity  estimates. 
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Thus,  we  suggest  that  the  remapping  error  plots  such 
as  that  shown  in  Fig.  11  can  be  used  to  segment  an  ice 
velocity  field  into  dark  regions  that  tend  to  be  large  rigidly 
rotating  plates  and  bright  regions  where  there  is  much 
shear  and  ice  deformation.  We  realize,  of  course,  that  this 
interpretation  is  based  on  a  single  observation.  Neverthe¬ 
less,  we  think  it  is  worth  pointing  out  since  it  may  well 
prove  to  be  a  useful  method  of  automatically  surveying 
sea-ice  velocity  fields  sensed  by  SAR  or  other  imaging 
sensors. 


VI.  Conclusions 

We  summarize  our  conclusions  based  on  the  research 
reported  above  as  follows: 

1)  For  the  sea-ice  images  considered  here  the  image 
pyramid  area  correlation  (IPAC)  algorithm  provides  a 
useful,  though  not  ideal,  means  of  automating  the  process 
of  deriving  sea-ice  velocity  information  from  SAR  im¬ 
ages.  Rejecting  low  confidence  velocity  estimates  and  re¬ 
placing  them  by  interpolated  values  improves  IPAC  re¬ 
sults  in  comparison  with  manual  results. 

2)  Feature  tracking  using  floe-lead  boundary  segments 
(BSFT)  produces  a  large  number  of  sea-ice  velocity  esti¬ 
mates  that  agree  well  with  IPAC  and  manual  results. 
However,  the  method  cannot  comprehensively  cover  the 
region  observed  by  the  SAR  because  trackable  floe-lead 
boundaries  are  not  present  everywhere.  Extension  to  pres¬ 
sure  ridge  features  is  indicated. 

3)  The  IPAC  and  BSFT  algorithms  provide  indepen¬ 
dent  estimates  of  sea-ice  velocity,  both  having  their  ad¬ 
vantages  and  disadvantages.  We  conclude  that  the  two  al¬ 
gorithms  complement  each  other  in  many  ways.  An 
algorithm  synthesized  from  both  techniques  is  indicated. 

4)  To  make  remotely  sensed  ice  velocity  field  data  use¬ 
ful  in  sea-ice  science  and  applications  a  measure  of  data 
quality  is  needed.  We  conclude  that  the  image  reconstruc¬ 
tion  and  subtraction  technique  provides  such  a  measure  of 
data  quality.  The  technique  may  also  be  useful  in  seg¬ 
menting  SAR  observations  of  sea  ice  into  rigid  plates  and 
regions  of  significant  shear  and  ice  deformation. 
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ABSTRACT 

Synthetic  aperture  radar  (SAR)  provides  an 
excellent  means  of  observing  the  movement  and 
distortion  of  sea  ice  over  large  temporal  and 
spatial  scales.  SAR  observations  are  not  affected 
by  clouds  or  lack  of  sunlight.  The  ERS-1  and 
JERS-1  satellites  will  carry  SAR's  over  polar 
regions.  A  key  component  in  utilization  of  these 
SAR  data  is  an  automated  scheme  for  extracting  the 
sea  ice  velocity  field  from  a  sequence  of  SAR 
images  of  the  same  geographical  region.  We 
describe  two  techniques  for  automated  sea  ice 
tracking,  image  pyramids  (hierarchical 
correlation)  and  feature  tracking.  Each  of  these 
techniques  is  applied  to  a  pair  of  SEASAT  SAR  sea 
ice  images.  The  results  compare  well  with  each 
other  and  manually  tracked  estimates  of  the  ice 
velocity  field.  Finally,  we  note  the  advantages 
and  di sadvantages  of  these  methods. 


1.  INTRODUCTION 

The  importance  of  sea  ice  remote  sensing  arises 
from  the  high  spatial  and  temporal  variability  of 
sea  ice,  its  wide  geographic  extent  and  its  impact 
on  climate,  arctic  off-shore  engineering, 
transportation  and  military  operations.  Sea  ice 
dramatically  changes  the  air-sea  energy  balance 
and  energy  transport.  Since  sea  ice  covers  some 
10X  of  the  Arctic  and  13%  of  the  Southern  Ocean, 
its  impact  upon  climate  is  strong. 

The  high  variability  of  sea  ice  coupled  with  its 
large  extent  makes  remote  sensing  techniques  a 
necessity  for  large-scale .studies.  The  ability  of 
synthetic  aperture  radar  to  provide  large,  high 
resolution  images  through  clouds  and  without 
sunlight  give  it  a  centralrole  in  sea  ice  remote 
sensing,  especially  ice  dynamics  (Ref.  1). 

Examples  of  SAR  sea  ice  images  are  given  in 
Fig.  1. 

The  process  of  retrieving  motion  information  from 
a  sequence  of  SAR  images  (of  the  same  geographic 
area)  revolves  around  two  principal  items,  namely 

1.  Identification  of  ice  features  to  use  as 
'tie  points'  for  ice  floe  tracking. 


2.  Finding  the  same  tie  points  in  subsequent 
images  of  the  sequence. 

This  process  has  been  accomplished  manually  using 
a  split-screen  video  display  and  an  experienced 
operator.  Ice  displacement  is  computed  auto¬ 
matically  once  the  locations  of  a  tie  point  are 
marked  in  a  pair  of  images  on  the  screen.  Results 
by  Thorndike  &  Rothrock  (in  Ref.  1)  using  this 
technique  are  given  in  Fig.  2  with  further  work  by 
Carsey  et  al.  and  Curlander,  Holt  and  Hussey  in 
Refs.  243 .  However,  this  process  is  too  slow  and 
expensive  for  real  time  and/or  extensive  use  with 
long-term  satellite  SAR's  such  as  ERS-1,  JERS-1 
and  RADARSAT.  This  paper  compares  two  approaches 
which  offer  good  prospects  for  automatically 
retrieving  sea  ice  displacement  information  from 
SAR  images.  The  development  of  an  automated 
technique  is  a  primary  requirement  for  successful 
use  of  satellite  or  aircraft  SAR  images  in  sea  ice 
remote  sensing. 


The  two  techniques  compared  here  are  called  the 
'image  pyramid*  and  ’boundary  segment’  sea  ice 
tracking  techniques.  The  most  advanced  of  the  two 
is  the  image  pyramid  or  hierarchical  correlation 
technique,  first  applied  to  sea  ice  by  Fily  and 
Rothrock  (Ref.  4).  Cross  correlation  of  blocks  of 
pixels  from  a  pair  of  cospatial  images  is  made 
computationally  efficient  by  first  building  from 
each  primary  high  resolution  image  a  sequence 
(pyramid)  of  derived  images  of  progressively 
coarser  resolution,  but  correspondingly  fewer  data 
points.  Block  cross  correlation  in  the  coarsest 
resolution  is  computationally  fast  and  builds  a 
crude  model  of  the  displacement  field  of  the  sea 
ice.  This* crude  model  is  used  to  guide  searches 
in  the  pyramid  of  progressively  higher  resolution 
images.  As  one  moves  d<xn  the  pyramid  to  the  " 
highest  resolution  image,  searches  are  guided  from 
previous  results  and  the  displacement  field  is 
continually  refined. 

The  second  technique  corresponds  more  closely  to 
the  manual  (visual)  technique.  Distinctive 
features  are  located  in  the  reference  (initial) 
image  and  then  located  again  in  the  test  (sub¬ 
sequent)  image.  This  is  a  pattern  recognition 
problem  where  techniques  are  dependent  upon  the 
image  content.  For  our  features  we  use  segments 
of  lead-floe  boundaries.  Thus  each  image  becomes 
a  set  of  boundary  segment  features.  The  set  of 
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Figure  1.  Two  SEASAT  SAR  images  of  the  same 

area  collected  over  the  Beaufort  Sea 
three  days  apart  on  orbits  1439  &  1482. 


features  in  the  reference  image  forms  the  tie 
points  which  are  then  located  in  the  test  image 
and  used  to  derive  the  ice  displacement  (velocity) 
field. 


2.  SEA  ICE  TRACKING  BY  IMAGE  PYRAMIDS 
2.1  Pyramid  Correlation 

The  problem  of  tracking  sea  ice  from  two  radar 
images  can  be  solved  by  finding  a  correspondence 
between  the  two  images  at  points  which  represent 
the  same  scene,  i.e.  the  same  sea  ice  image. 
Cross-correlation  of  blocks  of  pixels  in  the  images 
can  be  used  to  find  matching  scene  points  in  the 


images.  The  displacement  vector  between  matching 
points  is  taken  to  be  between  block  centers  with 
the  largest  correlation  coefficients.  This 
displacement  information,  together  with  the  tires 
of  acquisition  of  the  images,  can  be  used  to  derive 
velocity  information. 

Correlation  is  one  of  the  most  expensive  image 
processing  algorithms.  One  way  to  speed  up  the 
computation  is  to  apply  the  correlation  to  an  image 
pyramid  data  structure  (see  section  1).  Pyramids 
are  constructed  for  each  of  the  two  images  and  then 
the  following  algorithm  is  used: 

1.  For  each  non-overlapping  block  of  size 

(I  x  I)  in  the  coarsest  level  of  Pyramid  i, 
calculate  a  displacement  vector  by  finding 
the  global  max  of  the  correlation  with  the 
coarsest  level  of  Pyramid  2 

2.  For  each  finer  resolution  level  of  the 
pyramid  do  the  following: 

A.  Filter  the  displacement  vector  field  for 
the  previous  (coarser)  pyramid  level  in 
order  to  remove  gross  errors. 

B.  Interpolate  the  displacement  vector  field 
calculated  at  the  previous  pyramid  level 
to  provide  displacement  estimates  for 
each  (I  x  I)  non-overlapping  block  of  the 
current  level  of  Pyramid  1. 

C.  For  each  (I  x  I)  non-overlappi ng  block  of 
the  current  level  of  Pyramid  I  find  the 
local  max  of  the  correlation  with  the 
current  level  of  Pyramid  2  by  searching 
in  a  (M  x  M)  window  centered  at  the 
displacement  estimate  of  the  current 
level  of  Pyramid  2. 

3.  Optional:  Remove  displacement  vectors  with 

the  lowest  confidence  measure. 

The  algorithm  starts  at  the  top,  coarsest  pyramid 
level  and  proceeds  'top-down'  a  level  at  a  time. 

At  the  coarsest  level,  the  area  to  be  searched  in 
the  second  image  is  not  constrained.  It  is 
important  to  filter  gross  errors  in  the 
di splacement  vector  field  because  errors  in  the  top 
levels  of  the  pyramid  propagate  to  the  lower 
levels.  The  displacement  vector  field  is  then 
interpolated  between  pyramid  levels  since  an 
estimate  is  needed  for  those  blocks  for  which  there 
is  no  computed  displacement  at  the  coarser  level. 

At  each  level  in  the  pyramid,  the  estimated 
location  of  the  largest  correlation  found  at  the 
previous,  coarser  pyramid  level  constrains  the 
locations  searched  at  the  current  pyramid  level. 
These  location  constraints  decrease  the  total 
number  of  locations  that  must  be  searched. 


2.2  Experiments  with  SAR  ice  images 

The  two  SEASAT  Images  of  sea  ice  in  Fig.  1  were 
used  in  the  pyramid  correlation  experiments.  The 
images  were  first  averaged  in  10  x  10  blocks  to 
create  384  x  384  pixel  images.  The  averaged  images 
were  then  histogram  equalized.  Pyramids  with  three 
levels  were  created  from  both  images  by  using  a 
lowpass  filter  with  a  gaussian  transfer  function 
and  resampling.  The  resulting  levels  of  the 
pyramids  had  96  x  96,  192  x  192  and  384  x  384 
pixels. 
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Figure  2.  Sea  ice  velocity  field  obtained  by 
manual  tracking  of  Fig.  1  images. 


The  particular  implementation  of  the  above  pyramid 
correlation  algorithm  used  16  x  16  blocks  and 
32  x  32  search  windows  at  the  two  finer 
resolutions.  This  resulted  in  170  times  fewer 
calculations  than  direct  correlation.  The 
correlation  measure  used  in  all  the  experiments  was 
the  variance  normalized  cross-correl ation 
coefficient.  We  experimented  with  various  methods 
of  filtering  the  displacement  vector  field  in  order 
to  remove  gross  errors.  The  smoothing  filter  that 
provided  the  best  results  was  a  magnitude  median 
filter  where  each  displacement  vector  was  replaced 
by  the  vector  in  its  8  point  neighborhood  with  the 
median  magnitude.  This  filtering  method 
successfully  removed  displacements  with  errors  in 
magnitude,  but  did  not  perform  well  at  points  with 
errors  in  displacement  direction.  For  the  images 
used,  this  was  not  a  major  problem.  We  also  tried 
two  different  interpolation  methods,  a  nearest 
neighbor  interpolator  and  a  bilinear 
interpolator.  The  bilinear  interpolator  performed 
best  and  required  little  added  complexity. 


Figure  3.  Ice  displacement  vectors  generated 
by  the  image  pyramid  technique. 


The  experiment  using  lopass  pyramids  with  the 
median  filter  and  the  bilinear  i nterpol at i on  ga/e 
the  best  results.  Figure  3  shows  the  displacement 
vectors  at  the  384  x  384  resolution.  The  direction 
of  most  vectors  is  from  top  to  bottom.  Comparing 
these  displacement  vectors  with  the  manually 
generated  vectors  in  Figure  2,  it  is  seen  that  the 
method  performs  very  well.  The  problem  areas  at 
the  edges  of  the  figure  correspond  to  ice  features 
flowing  into  or  out  of  the  finite  area  covered  by 
the  images.  Many  errors  are  removed  by  choosing 
only  those  displacement  vectors  with  high 
correlation  coefficients.  Fig.  4  shows  vectors 
where  the  correlation  coefficient,  C  >  0.4. 


Figure  4.  Ice  displacement  vectors  by  image  pyramid 
with  low  confidence  vectors  removed. 


The  errors  in  the  interior  of  the  image  field  cause 
more  concern.  From  a  qualitative  evaluation  it 
seemed  that  some  of  these  errors  occurred  near 
leads,  where  the  shapes  of  ice  features  change  so 
much  that  the  correlation  maximum  shifts  to  some 
erroneous  location.  One  of  the  errors  occurred  in 
the  middle  of  an  ice  field  which  had  very  little 
structure,  and  thus  resulted  in  a  broad  correlation 
peak.  We  assume  that  this  peak  was  swamped  out  by 
noise  or  systematic  distortions  such  as  rotation. 

2.3  Pi scussion 

The  properties  of  the  pyramid  correlation  algorithm 
make  it  well  suited  to  automatic  tracking  of  ice 
floes.  The  search  originally  takes  place  at  a 
coarse,  global  level.  This  is  analogous  to  people 
scanning  the  images  to  find  a  general  mafch.  The 
search  at  the  finer  levels  is  analogous  to  people 
focusing  their  attention  to  specific  areas  for  a 
more  accurate  match.  The  algorithm  provides  noise 
immunity  by  starting  the  search  at  the  coarse 
level,  where  averaging  has  reduced  susceptibility 
to  noise.  The  constraints  at  the  finer  levels 
provide  noise  immunity  by  not  allowing  a  search  at 
all  locations.  Limiting  the  search  locations 
lowers  the  probability  of  detection  of  a  false  peak 
due  to  noise.  Finally,  the  algorithm  is  much  more 
efficient  than  direct  correlation. 

The  correlation  measure  currently  used  cannot 
detect  features  which  undergo  large  rotations. 

This  may  be  alleviated  by  using  a  rotationally 
invariant  measure  currently  being  developed. 
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3.  SEA  ICE  TRACKING  USING  FLOE-LEAD 
BOUNDARY  FEATURES 

3.1  Introduction  to  the  approach 


The  second  sea  ice  tracking  algorithm  is  a  pattern 
recognition  approach.  Here  we  have  chosen  floe- 
lead  bounary  segments  as  our  physical  feature.  The 
movement  of  these  features  from  one  image  to  the 
next  provides  the  information  from  which  the  sea 
ice  velocity  field  is  estimated.  Feature  patterns 
must  be  very  distinctive,  if  not  unique,  objects. 
They  are  character! zed  by  information  which  is 
Independent  of  the  object  location  and  orientation, 
though  this  latter  information  is,  of  course, 
attached  to  them. 

Thus  each  image  becomes  a  set  of  features  described 
by  parameters  independent  of  the  feature  location 
and  orientation.  The  collection  of  several  hundred 
features,  each  described  by  several  hundred 
numbers,  reduces  the  information  needed  to 
characterize  a  100  X  100  km  image.  The  search 
space  is  correspondingly  reduced. 


Under  many  circumstances  ice  and  water  in  SAP  irage 
can  be  distinguished  by  a  simple  pixel-intensity 
threshold  (Ref.  A).  We  use  such  a  threshold  to 
reduce  the  SAR  image  to  binary  pixels,  1  for  floe 
ice  and  0  for  lead  water  (including  grease  ice). 
Floe-lead  boundary  edges  are  identified  by  simply 
finding  water  pixels  adjacent  to  ice  pixels. 

* 

Next  the  water  edge  pixels  identified  above  are 
classified  into  eight  types  according  to  their 
surroundings  to  determine  the  orientation  of  the 
floe-lead  boundary  edge.  The  four  nearest  neighbor 
pixels  are  considered  to  find  if  there  is  water 
above,  below,  to  right  or  to  left  of  a  given 
boundary  pixel.  This  scheme  Is  Illustrated  below. 

O 


O 
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Once  a  feature  pattern  is  recognized  in  a  pair  of 
images,  sea  ice  velocity  information  can  be 
derived.  Since  we  know  something  of  a  feature's 
physical  characteristics,  we  can  derive  not  just 
one  but  many  velocity  vectors,  knowing  the 
displacement  and  rotation  of  a  particular 
feature.  That  is,  we  know  the  movement  of  every 
pixel  in  the  feature.  The  results  reported  here 
are  intended  as  a  progress  report  on  this  pattern 
recognition  approach  to  sea  ice  tracking. 
Refinement  of  the  technique  is  in  progress. 


3,2  Pattern  recognition  tracking  with  floe  lead 
boundary  segments 

We  first  present  the  boundary  segment  tracking 
algorithm  in  outline  form. 

1.  SAR  ice  field  images  averaged  to  obtain 
250  x  250  m  pixels. 

2.  SAR  image  transformed  into  binary  (ice  or 
lead)  image  by  thresholding. 

3.  Edge  detection  algorithm  transforms  image 
into  pattern  of  floe-lead  boundaries. 

4.  Boundary  pixels  classified  according  to 
surroundings. 

5.  Boundary  segments  assembled  from  classified 
pixels  such  that  segment  corresponds  to  only 
one  side  of  a  lead. 

6.  Segments  characterized  by  deviation  from 
straight  line  fit. 

7.  Most  distinctive  features  in  reference  image 
are  searched  for' in  test  image  to  give 
initial  estimate  of  ice  movement. 

8.  Initial  estimate  of  ice  velocity  field  guides 
search  for  less  distinct  features. 

9.  Each  segment  successfully  found  in  test  image 
yields  many  ice  velocity  vectors. 


Figure  5.  Diagram  illustrating  edge  pixel 
classification  scheme. 


The  classified  pixels  are  then  assembled  into 
boundary  segments.  The  assembly  scheme  used  here 
begins  by  finding  the  dominant  direction.  For  the 
orbit  1439  image  the  statistics  for  edge  pixel 
direction  clearly  indicated  that  the  dominant  lead 
direction  is  north-south,  and  similarly  for  the 
1482  image.  Consequently  we  divided  our  boundary 
segments  into  two  groups  corresponding  to  the  left 
and  right  edges  of  the  leads,  to  wit: 

Left  edges  (types  2,3,4  and  associated  1  &  5) 

Right  edges  (types  6,7,8  and  associated  1  &  5) 

Our  aim  here  is  to  match  left  edges  in  image  1439 
with  left  edges  in  image  1482  and,  independently, 
right  edges  in  the  image  pair.  This  allows  ice 
floes  along  the  left  side  of  a  lead  to  slip  with 
respect  to  ice  floes  along  the  right  side.  Thus 
shear  motion,  expected  across  leads,  is  built  into 
the  algorithm. 

Boundary  segment  features  are  characterized  by 
their  deviation  from  a  straight  line.  First,  a 
"least  squares"  criterion  straight  line  is  fit  to  a 
given  boundary  segment.  The  distance  from  this 
straight  line  to  the  segment  is  calculated  at 
regular  intervals  along  the  straight  line.  This 
set  of  deviations  characterizes  the  shape  of  a 
boundary  segment  as  illustrated  below.  Note  that 
this  method  of  characterizing  a  boundary  segment 
feature  is  invariant  under  rotation  of  the  feature. 


3.3  Description  of  the  algorithm 

Since  SAR  pixels  cannot  be  geographically  located 
better  than  about  200  m  and  this  accuracy  is 
adequate  for  scientific  and  practical  applications, 
we  average  10  x  10  SAR  pixel  blocks  to  obtain  an 
image  with  about  250  m  resolution.  This  reduces 
computation  time  and  SAR  speckle  noise. 


Figure  6.  Edge  segment  fitted  by  straight  line 

and  characteri zed  by  deviation  from  line. 


Recognition  of  matching  feature  patterns  is  done 
simply  by  comparing  the  set  of  deviation  measure¬ 
ments  for  two  candidate  features.  A  variance 
normalized,  cross  correlation  coefficient  (C)  was 
used  to  evaluate  how  closely  two  candidate  features 
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matched.  Often  segments  which  we  want  to  compare 
are  of  different  lengths  because  of  partial  break¬ 
up  of  a  feature  between  two  successive  SAR 
Images.  So  in  our  search  we  include  features  which 
are  of  differing  lengths.  The  shorter  feature  is 
slid  along  the  longer  feature  and  C  is  calculated 
for  all  placements  of  the  smaller  feature  with 
respect  to  the  larger.  This  allows  for  large 
features  to  break  Into  parts  and  still  be  found. 

We  have  at  this  point  two  images  which  each  contain 
two  sets  of  boundary  segment  features  -  left  and 
right  edges.  The  search  procedure  must  select  a 
feature  from  the  reference  image  (first  in  time), 
and  try  to  match  it  to  a  feature  in  the  test  image 
(later  in  time).  A  successful  match  reveals  ice 
velocity  vectors  all  along  the  matching  segments. 

The  search  and  matching  procedure  begins  by 
considering  only  the  larger  edge  features.  These 
features  are  more  likely  to  be  correctly  matched 
and  because  of  the  large  size  reveal  the  ice 
velocity  field,  ?(x,y)  ,  over  a  large  region.  For 
a  left  edge  feature  with  n  pixels  we  search  over 
left  edge  features  having  lengths  from  n/2  to 
2n  pixels.  The  search  is  restricted  in  area  to 
displacements  less  than  60  km  -  corresponding  to  a 
maximum  ice  movement  of  20  km  per  day  over  the 
three  days  between  the  1439  and  1482  images. 

Correct  matches  for  both  the  sets  of  left  and  right 
edges  give  two  sets  of  ice  displacement  vectors. 

To  eliminate  erroneous  matches  we  disregard 
velocity  vectors  which  do  not  come  reasonably  close 
to  their  neighbors  in  both  magnitude  and 
direction.  This  set  of  matches  between  large 
features  yields  an  initial  estimate  of  the  ice 
velocity  field  which  can  then  be  used  to  guide  the 
search  for  matches  between  shorter  features. 
Successful  identification  of  a  feature  in  both  the 
reference  and  test  images  yields  much  ice  velocity 
information.  Since  a  feature  is  generally  some 
tens  to  even  hundreds  of  pixels  in  length,  many  ice 
velocity  vectors  can  be  derived  by  displacement  of 
each  pixel  in  the  shorter  of  the  matching  pair  of 
features. 


3.4  Application  to  SEASAT  SAR  images  in  the 
Beaufort  Sea 

We  now  apply  the  boundary  segaent  algorithm  to  the 
images  of  Fig.  1.  Since  the  algorithm  is  still 
being  developed,  the  results  presented  here  are 
only  for  a  first  pass  through  the  image  data. 
Boundary  edges  constructed  from  the  images  of 
Fig.  1  varied  In  length  from  4  to  163  pixels.  For 
this  first  pass  we  considered  only  the  longer 
edges.  Thus  in  the  reference  image  (orbit  1439) 
edges  of  30  or  more  pixels  were  used.  In  the  test 
image  (orbit  1482)  matching  segments  of  length 
(n/2)  »  (30/2)  *  15  pixels  and  above  were 
considered.  The  resulting  sets  of  boundary 
segments  for  the  orbit  1439  and  1482  images  are 
shown  below. 

Each  1439  edge  was  compared  to  the  set  of  1482 
edges  seeking  the  best  match  by  maximizing  C.  The 
search  was  restricted  to  edges  within  60  km  of  the 
1439  edge  location  since  ice  movement  of  over  20 
km/day  is  unlikely.  After  the  displacement  vectors 
(between  segment  centers)  were  calculated  for 
matching  left  edges  in  the  two  images,  this  process 
was  repeated  independently  for  right  edges.  Only 
then  were  results  compared. 


Figure  7.  Long  boundary  segments  constructed  from 
orbit  1439  (left)  and  orbit  1482  (right)  SAR 
image  data. 

On  this  first  pass  through  the  data  we  want  to 
insure  that  only  correct  displacement  vectors  are 
selected.  Thus  we  demanded  rather  strict  agreement 
between  neighboring  displacement  vectors.  We  thus 
assumed  that,  for  a  small  area  of  the  image,  the 
correct  left  and  right  (edge)  displacement  vectors 
should  not  vary  in  length  by  more  than  30*  or  have 
an  angular  deviation  greater  than  10  degrees.  Only 
those  left  edge  displacement  vectors  for  which  a 
similar  right  edge  displacement  vector  could  be 
found  within  15  km  were  recorded.  Although  this 
constraint  eliminated  a  few  correct  displacements, 
it  also  eliminated  many  erroneous  displacements 
arising  from  correlations  between  nearly  straight 
edges. 

Fig.  8  below  shows  the  displacement  vectors  for  the 
centers  of  the  matched  boundary  segments.  Each  of 
the  vectors  represents  a  family  of  displacement 
vectors  resulting  from  the  known  displacement  of 
each  pixel  in  the  segment  pair.  A  more  dense  and 
extended  field  of  ice  displacement  vectors  is  thus 
Summarized  by  these  few  vectors. 


Figure  8.  Ice  displacement  (velocity)  vectors 

from  matching  of  long  boundary  segments  in 
images  1439  and  1482.  Solid  lines  are  left 
edges  and  dotted  lines  right  edges. 


These  displacement  vectors  agree  well  with  the 
manually  generated  vectors  of  Fig.  2  and  the  image 
pyramid  generated  vectors  of  Figs.  3  and  4.  From 
these  initial  results  it  is  clear  that  successive 
passes  will  have  to  be  made  through  the  data.  The 
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displacement  vectors  calculated  on  this  first  pass 
can  be  used  to  narrow  the  search  area  on  subsequent 
passes  with  the  goal  of  finally  calculating  a 
displacement  vector  for  a  majority  of  the  edges. 
Other  distinctive  sea  ice  features  such  as  pressure 
ridges  could  also  be  used  in  feature  tracking  just 
as  the  boundary  segments  are  here. 


4.  COMPARISON  OF  METHODS 

Given  our  results,  we  can  qualitatively  address 
some  of  the,  advantages  and  di sadvantages  of  the 
above  techniques.  In  particular,  we  can  compare 
the  methods  using  an  outline  suggested  by  Kashef 
and  Sawchuck  (Ref.  5). 

A  feature-matching  method,  such  as  'boundary 
segment/  approaches  map  (image)  data  differently 
from  a  correlation  method,  such  as  'image 
pyramid/  Whereas  correlation  methods  rely  on  the 
intensity  values  of  map  data,  feature-matching 
methods  rely  on  detected  edges,  boundaries,  and 
vertices  of  map  data.  The  advantage  of  using 
'boundary  segment'  in  tracking  sea  ice  lies  in 
reducing  an  image  of  millions  of  pixels  to  one  of 
thousands  of  pixels.  Only  relevant  data  applies  to 
tracking.  The  potential  cost-savings  for 
computational  speed  exist,  depending  on  how 
'relevant*  the  reduced  data  is.  In  particular,  one 
should  only  pick  those  edges,  boundaries,  and 
vertices  common  to  both  images  for  tracking. 

Picking  those  edges,  boundaries,  and  vertices 
common  to  both  images  remains  difficult.  Thus,  the 
advantage  for  'boundary  segment'  also  harbors  a 
disadvantage.  Boundaries  may  fragment,  rotate  or 
disappear.  Furthermore,  detectable  boundaries  may 
depict  more  than  one  floe,  thereby  compounding  the 
tracking  problem.  Hence,  the  disadvantage  to 
'boundary  segment*  is  the  difficulty  in  identifying 
those  edges,  boundaries,  and  vertices  conroon  to 
both  images. 

Correlation  methods  surmount  this  difficulty  by 
assuming  little  or  nothing  about  image  content. 
Traditionally  computation-intensive,  correlation 
methods  become  attractive  with  some  of  the  newer 
algorithms  such  as  pyramid  correlation.  As 
discussed,  'image  pyramid*  techniques  can  track 
sea-ice  floes  in  a  manner  more  efficient  than 
direct  correlation.  The  noise  immunity  provided  by 
'image  pyramid'  techniques  comes  as  an  added 
bonus.  Therefore,  the  advantage  of  using  'image 
pyramid'  techniques  stems  threefold:  'image 
pyramids'  assume  little  about  image  content;  'image 
pyramids'  can  track  global  displacements  well; 
'image  pyramids'  have  inherent  noise  immunity. 

Oespite  the  advantages  of  using  ’image  pyramid' 
techniques,  'image  pyramid*  techniques  like  the  one 
currently  used  for  this  research  cannot  detect 
small-scale  pattern  changes  which  involve  a  large 
rotation.  This  problem  is  significant  given  the 
rotational,  translational,  nature  of  sea  ice. 

To  address  this  drawback,  one  can  opt  for  two 
schemes.  One  scheme  involves  substituting  a 
rotation-invariant  measure  in-place  of  standard 
correlation.  Another  scheme  synthesizes  a  hybrid 
method  from  feature-tracking  and  correlation.  For 
example,  the  feature-matching  method  presented  here 
is  conditioned  to  small-scale  rotation  changes.  A 
synthesis  of  feature-tracking  and  correlation 
methods  seems  likely  in  the  future. 


5.  SUMMARY  AND  CONCLUSIONS 

Below  *e  summarize  major  results  and  conclusions. 
Significant  further  research  is  needed  to  produce 
an  efficient  and  reliable  algorithm  for 
automatically  extracting  the  sea  ice  velocity  field 
from  SAR  image  data.  ' 

1.  Estimates  of  the  sea  ice  velocity  field 
are  successfully  -derived  from  SAR  images  using  the 
image  pyramid  (hierarchical  correlation)  technique. 

2.  Further  refinement  of  the  image  pyramid 
technique  should  include  removal  of  a  small  number 
of  erroneous  velocity  vectors  using  constraints 
dependent  on  correlation  coefficient,  near  neighbor 
velocities,  velocity  gradients,  etc. 

3.  Estimates  of  sea  ice  velocity  are 
successfully  obtained  over  limited  areas  using  the 
boundary  segment  tracking  technique. 

4.  Further  refinement  of  the  feature  tracking 
technique  should  include  use  of  pressure  ridge 
features,  constrained  search  and  relaxation 
techniques. 
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ABSTRACT 

Synthetic  aperture  radar  (SAR)  provides  an  excellent 
means  of  observing  the  movement  and  distortion  of  sea 
ice  over  large  temporal  and  spatial  scales.  SAR 
observations  are  not  affected  by  clouds  or  lack  of  sunlight. 
The  European  Space  Agency's  ERS-1  satellite  as  well  as 
others  planned  to  follow  will  carry  SAR's  over  the  polar 
regions  beginning  in  1989.  A  key  component  in  utilization 
of  these  SAR  data  is  an  automated  scheme  for  extracting 
the  sea  ice  velocity  field  from  a  sequence  of  SAR  images 
of  the  same  geographical  region.  Several  schemes  are 
considered.  The  image  pyramid  area  correlation  (IPAC) 
method,  also  known  as  hierarchical  correlation,  has 
successfully  performed  automated  sea  ice  tracking,  but  is 
vulnerable  to  uncertainties  if  the  ice  rotates  more  than 
about  1 0°  to  1 5°  between  SAR  observations.  Since  ice 
motion  is  not  uniform,  window  size  is  important  during 
cross  correlation  operations.  Selection  of  optimum 
window  size  depends  on  requirements  regarding 
accuracy  of  location  and  suppression  of  spurious 
sidelobes.  Rotation-invariant  methods  have  the  potential 
to  succesfully  track  isolated  floes  in  the  marginal  ice  zone 
and  may  also  be  helpful  in  cases  of  rotating  pack  ice. 

Hu's  invariant  moments  are  considered  as  a  possible 
candidate  for  rotation-invariant  tracking.  Feature  tracking 
is  inherently  robust  with  respect  to  tracking  rotating  sea 
ice.  but  has  some  limitations  when  floe-lead  boundaries 
are  used  as  features.  The  use  of  pressure  ridge  features 
is  indicated.  It  appears  that  an  automated  sea  ice  tracking 
algorithm  will  need  a  variety  of  techniques  to  operate 
successfully  over  a  wide  range  of  sea  ice  conditions. 

Keywords:  Sea  ice,  synthetic  aperture  radar  (SAR), 
automated  tracking,  rotation,  image  pyramid,  feature 
tracking,  invariant  moments 


I.  INTRODUCTION 

Sea  ice  remote  sensing  is  important  because  of  the 
high  spatial  and  temporal  variability  of  sea  ice,  its  wide 
geographic  extent  and  its  impact  on  weather,  climate, 
arctic  off-shore  engineering,  transportation  and  military 


operations.  Sea  ice  dramatically  changes  the  air-sea 
energy  balance  and  energy  transport.  Since  sea  ice 
covers  significant  and  variable  areas  of  the  Arctic  and 
Southern  Oceans,  it  has  a  significant  impact  on  both  polar 
and  global  weather  and  climate. 

The  high  variability  of  sea  ice  coupled  with  its  large 
extent  makes  remote  sensing  techniques  a  necessity  for 
large-scale  ice  science  as  well  as  for  timely  practical 
applications,  such  as  deployment  of  ocean  drilling 
platforms.  Since  polar  regions  are  often  cloud  covered 
and/or  without  sunlight  for  long  periods,  the  ability  of 
synthetic  aperture  radar  (SAR)  to  provide  large  scale,  high 
resolution  images  through  clouds  and  without  sunlight 
gives  it  a  central  role  in  sea  ice  remote  sensing.  SAR 
provides  good  information  on  sea  ice  extent,  movement 
and  deformation,  internal  geometry  (structure  of  floes  and 
leads)  and  surface  roughness  [1].  Some  useful 
information  on  sea  ice  types,  thickness,  and  physical 
properties,  can  also  be  obtained  from  SAR  images. 
Examples  of  SAR  images  of  sea  ice  are  given  in  Fig.  1 
below  and  by  Weller  et  al.  [1]. 

The  process  of  retrieving  ice  motion  information  from 
a  sequence  of  SAR  images  (of  the  same  geographic  area) 
revolves  around  two  principal  items,  namely: 

1 .  Identifying  regions  in  a  sea  ice  field  to  use 
as  “tie  points'  for  ice  floe  tracking. 

2.  Finding  the  same  tie  points  in  subsequent 
images  of  the  sequence. 

Once  ice  displacement  is  known,  ice  velocity  can  easily 
be  calculated  since  the  elapsed  time  between- the  SAR 
images  is  known.  In  this  paper  we  shall  consider  the  ice 
velocity  and  displacement  fields  as  equivalent  since  they 
are  simply  related  by  a  scale  factor. 

The  tracking  process  has  been  accomplished 
manually  using  a  split-screen  video  display  and  an 
experienced  operator.  Ice  velocity  is  computed 
automatically  once  the  locations  of  tie  points  are  marked 
in  a  pair  of  images  on  the  screen.  Thorndike  &  Rothrock 
[1]  have  produced  interesting  results  using  this  technique; 
further  work  has  been  done  by  Carsey  et  al.  (2)  and  by 
Curlander,  Holt  and  Hussey  [3].  However,  this  process  is 
too  slow  and  expensive  for  real  time  and/or  extensive 
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Fig- 1.  A  pair  of  co-spatial  SEASAT  SAR  images  of  sea  ice  in  the 
Beaufort  Sea  from  orbits  1396A  (top)  and  14390  (bottom)  ,  collected 
on  October  2  and  Octobers.  1978.  Banks  Island  lies  along  the  right 
edge  of  the  images. 

analysis  of  images  provided  by  long-term  satellite  SARs. 
Two  and  possibly  three  such  satellite  SARs  are  scheduled 
to  be  in  orbit  by  the  early  1990s,  namely  the  European 
ERS-1 ,  the  Japanese  ERS-1  and  the  joint  Canadian  / 
U.S.A.  RADARSAT.  Rea!  aperture  radar  images  at  low 
resolution  are  apparently  also  available  from  the  Soviet 
Kosmos  1 500  series  satellites.  Hence,  the  development 
of  automated  techniques  is  a  primary  requirement  for 
successful  use  of  satellite  or  aircraft  SAR  images  in  sea 
ice  remote  sensing. 

Automated  ice  tracking  algorithms  have  been 
developed  which  work  reasonably  well  in  regions  where 
there  is  little  rotation  and  distortion,  e.g.  in  portions  of 


central  pack  ice  [4  and  5].  However  in  regions  where 
there  is  significant  ice  rotation,  greater  than  1 0'  to  1 5°,  the 
existing  image  pyramid  area  correlation,  IPAC,  (also 
called  hierarchical  correlation)  methods  are  not  expected 
to  work  well  because  the  cross  correlation  comparisons 
are  made  between  image  blocks  which  have  been 
translated  in  position,  but  not  rotated.  For  a  much  used 
image  pair  (SEASAT  orbitsl  439  and  1482)  upon  which 
the  image  pyramid  algorithm  works  well,  the  maximum  ice 
rotation  is  some  5°.  Recent  work  by  Vesecky  et  al.  [6]  has 
shown  that  ice  motion  containing  shear  and  distortion  can 
also  cause  errors  when  using  image  pyramid  methods. 

In  this  paper  we  discuss  the  impact  of  ice  rotation  and 
distortion  on  sea  ice  tracking  algorithms  and  report 
several  new  results  concerning  the  performance  of 
existing  algorithms  and  methods  for  improving  tracking 
performance  when  significant  ice  rotation  and  distortion 
are  present.  In  section  2,  we  discuss  the  effects  of  rotation 
and  distortion  on  the  IPAC  algorithm.  In  this  section  we 
report  measurements  of  block  cross-correlation  function 
variations  with  respect  to  rotation  angle  and  shear  angle. 
The  effects  of  window  size  on  cross  correlation  between 
image  pairs  are  also  discussed.  In  section  3,  we  consider 
the  use  of  rotational  invariant  measures  to  overcome  the 
vulnerability  of  the  IPAC  method  to  ice  rotation.  This 
method  is  applied  to  ice  floes  which  have  rotated  during 
the  time  interval  between  collection  of  the  members  of  an 
image  pair.  A  case  study  using  SAR  images  from 
SEASAT  orbits  1438  and  1481  is  examined.  In  section  4, 
we  discuss  briefly  the  advantages,  and  some  limitations, 
of  using  our  feature  tracking  algorithm  [5,6]  to  track 
rotating  ice  fields.  We  summarize  our  conclusions  in 
section  5. 


2.  SENSITIVITY  OF  BLOCK  CORRELATION  TO 
SEA  ICE  ROTATION  AND  DISTORTION 

Block  correlation  for  matching  in  a  pyramid  data 
structure  has  been  applied  to  SEASAT  sea  ice  scenes 
[4,5,6].  We  have  sectioned  pairs  of  SAR  images  of  the 
same  geographic  location  into  non-overlapping  image 
pairs.  The  IPAC  algorithm  has  been  successfully  applied 
in  these  pack  ice  regions.  It  has  also  performed 
successfully  on  the  pair  1439  and  1482  taken  as  a  whole. 

Pure  translational  motion  is  the  strict  assumption  one 
makes  when  the  correlation  coefficient  is  used  as  a 
similarity  measure  between  images.  Since  pack  ice 
regions  are  observed  to  move  as  large  blocks, 
approximating  rigid  translation,  we  believe  that  IPAC  will 
sucessfully  track  these  regions.  Here  we  assess  the 
performance  of  the  algorithm  in  the  presence  of  rotation 
and  non-rigid  motions.  These  types  of  motion  are 
expected  in  the  marginal  ice  zone. 

The  study  of  the  degradation  of  the  correlation 
coefficient  due  to  geometric  distortions  is  important  for  two 
reasons.  First,  it  can  be  used  to  determine  the 
applicability  of  correlation  when  the  motion  is  not  strictly 
translation.  Second,  it  can  be  used  to  set  suitable 
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Fig  2  Perspective  plots  of  the  two  dimensional  cross-correlation  function  between  matching  SAR 
sea  ice  scenes  from  SEASAR  orbits  1439  and  1482.  Image  resolution  was  averaged  to  250  meters.  The 
window  size  is  increases  in  steps  of  a  factor  2  moving  from  left  to  right:  i.e.  window  sizes  are  5,  9. 17,  33,  65, 
and  129  pixels  to  a  side. 


correlation  window  sizes. 

It  is  known  that  the  linear  correlation  coefficient  is 
affected  by  noise  and  by  geometric  distortions  (7J.  The 
noise  considerations  apply  to  all  applications  of 
correlation,  including  registration,  where  it  is  most 
commonly  used.  Our  application,  however,  is  to  motion 
detection  and  the  geometric  distortions  that  occur  are 
basically  the  result  of  nonuniform  motion  and  are  an 
expected  occurrence.  In  this  section  we  present  results 
related  mainly  to  geometric  distortions. 

The  analysis  of  the  degradation  of  correlation 
coefficient  due  to  geometric  distortions  is  complex  and, 
unfortunately,  data  dependent.  For  this  reason,  we  have 
conducted  an  experimental  appraisal  of  distortion  effects 
on  real  ice  data. 


2.1  Scale  of  the  Motion 

If  the  motion  between  the  two  images  is  pure 
translation,  then  one  can  in  principle  find  the  motion  by 
cross-correlating  the  two  images.  In  statistical  terms,  the 


problem  is  to  find  an  estimator  for  the  translation.  In  order 
to  find  a  minimum  variance  estimator,  given  noisy  data, 
one  would  use  all  the  data  one  has  and  this  results  in 
chosing  the  largest  window  possible. 

If  the  motion  between  the  two  images  is  not  a  pure 
translation,  then  there  is  a  tradeoff  involved  in  picking  a 
window  size  for  correlation.  Larger  windows  are  desired 
for  accuracy  in  the  estimators,  but  are  undesirable  for 
scales  of  motion  smaller  than  the  window  size.  Smaller 
windows,  on  the  other  hand,  are  desirable  if  the  motion  is 
only  locally  uniform  but  are  undesirable  for  accuracy  of 
the  estimators,  especially  in  the  presence  of  noise. 

We  show  in  Figure  2  a  sequence  of  perspective  plots 
of  the  cross  correlation  of  a  matching  section  of  1439  and 
1 482  which  were  averaged  to  250  meter  pixels.  From  left 
to  right  in  the  figure,  the  window  size  is  increased  by  a 
factor  of  about  two  at  each  step.  At  small  window  sizes, 
the  cross-correlation  is  evidently  noisy  and  may  result  in 
misses  if  a  large  enough  region  is  searched  for  a  match. 
We  find  the  window  size  1 7x1 7  to  have  a  distinct  peak 
and  it  is,  therefore,  suitable  for  use.  The  larger  window 
sizes  show  a  decrease  in  correlation  coefficient.  We 
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observe  nonuniform  motion  in  this  region  of  the  images 
and  we  attribute  the  decrease  in  the  correlation  coefficient 
to  this  non-uniformity.  Although  the  peak  for  these  larger 
window  sizes  may  be  easily  extracted,  it  would  not 
represent  the  motion  accurately. 

V.  E.  Borodachev  [8]  has  studied  aerial  photographs 
and  has  found  the  average  size  of  unfractured  ice  to  be  10 
kilometers.  This  is  about  the  size  of  the  1 6x1 6  windows 
we  currently  use  and  further  justifies  the  window  size  if  we 
assume  that  one  or  two  ice  blocks  move  uniformly. 


2.2  Rotation  and  Shear 

We  have  studied  the  effect  of  rotation  and  shear  by 
simulating  distortions  on  the  computer.  We  are  able  to 
remap  an  image  by  specifying  any  linear  transformation 
between  the  coordinates  of  the  original  image  and  the 
output  image  and  then  using  bicubic  interpolation  to 
resample  the  image.  In  particular,  we  can  simulate 
rotation  and  pure  shear  using  the  following  matrices: 

Rotation 

(cos  0  -sin  fl\ 
sin  0  cos  0  ) 

Shear 

1  tan  0\ 

0  1  ) 


Our  experimental  study  consisted  of  extracting  twenty 
non-overlapping  16x16  subimages  at  125  meter  pixels  of 
1439  and  1396  and  applying  both  transformations  to  each 
subimage  at  increments  of  one  degree.  We  then  find  the 
cross-correlation  between  the  original  and  the  distorted 
image  for  a  window  size  of  16x16  pixels,  which  is  used  in 
our  current  implementation  of  the  IPAC  algorithm  (5,6). 
Figure  3  shows  a  plot  of  the  average  correlation 
coefficient  for  the  twenty  subimages  as  a  function  of  angle 
of  rotation  or  angle  of  shear,  as  defined  by  the  matrices 
above. 

If  the  correlation  coefficient  of  the  correct  match  is 
low,  then  the  probability  of  a  mismatch  increases.  In  one 
case,  when  we  searched  a  local  area  around  the  correct 
match,  we  found  a  higher  correlation  coefficient  than  the 
correct  one  when  the  rotation  angle  was  1 6  degrees.  The 
maximum  correlation  coefficient  was  .39  and  the  distance 
to  the  correct  match  was  about  3  pixels.  We  believe  that 
rotations  of  ten  degrees  will  not  hamper  detection  by 
block  correlation. 

One  can  intuitively  interpret  the  fact  that  the  shear 
curve  lies  above  the  rotation  curved  by  considering  the 
overlap  of  pixels  for  a  given  shear  or  rotation.  The  pixels 
in  the  middle  row  of  the  image  are  fixed  for  a  shear 
whereas  only  the  middle  pixel  is  fixed  for  a  rotation. 

A  further  study  determines  the  effect  of  resolution  on 
rotation  effects.  Figure  4  shows  the  effects  of  rotation  on 
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Fig.  3.  Averaged  two-dimensional  cross-correlation  function  with 
respect  to  angle  of  rotation  and  angle  of  shear.  Twenty  SAR  sea  ice 
scenes,  extracted  from  SEASAT  orbits  1396  and  1439,  were  used  in 
the  averaging  process. 

images  depicting  the  same  scene  at  various  levels  of 
resolution,  corresponding  to  the  pyramid  levels  used  in 
the  current  IPAC  algorithm.  We  note  that  the  coarsest 
level  of  resolution  is  more  immune  to  rotation  effects  than 
the  equivalent  window  would  be  at  the  finest  level.  This  is 
an  advantage  of  pyramid  correlation  in  areas  of  rigid 
rotation. 

In  summary,  correlation  is  strictly  valid  only  for 
translational  motion.  Studies  of  the  actual  ice  data,  as 
represented  by  Figure  2,  have  shown  that  the  ice  motion 
is  not  translational,  but  has  rotational  and  nonuniform 
translational  components.  We  have  shown,  however,  that 
correlation  is  still  practical  under  the  various  distortions 
studied.  Studies  of  simulations  of  rotations  and  shears,  as 
depicted  in  Figure  3,  have  shown  that  correlation  is  useful, 
up  to  moderate  angles  of  distortion.  In  particular,  we 
believe  that  rotations  of  about  10°  can  be  dealt  with 
adequately  with  correlation.  Finally,  Figure  4  shows  the 
effect  of  correlation  in  a  pyramid.  As  expected,  smaller 
windows,  at  the  coarse  level  are  more  immune  to  rotation 
that  their  equivalent  windows  could  be  at  the  fine 
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Fig.  4.  Two-dimensional  cross-correlation  function  with  respect  to 
angle  of  rotation  for  window  sizes  of  16, 32,  and  64  pixels  to  a  side  all 
depicting  the  same  area  of  16x16  kilometers.  Twenty  SAR  sea  ice 
scenes,  extracted  from  SEASAT  orbits  1396  and  1439,  were  used  in 
the  averaging  process. 
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resolution  level. 


3.  APPLICATION  OF  INVARIANT  MOMENTS 

FOR  TRACKING  SEA  ICE  WHICH  ROTATES 

3.1  Invariant  Moments:  Background 

Invariant  moments  help  solve  the  problem  of 
identifying  objects  in  a  scene  that  changes  because  the 
objects  rotate  from  one  image  of  the  scene  to  the  next.  As 
mentioned  in  the  previous  section,  IPAC  will  handle 
rotations  in  sea  ice  features  up  to  about  10°  to  15°. 

Beyond  that,  rotations  will  start  to  affect  the  ability  of  the 
algorithm  to  correctly  track  features.  Rotations  of  that 
magnitude  or  greater  commonly  occur  in  places  like  the 
marginal  ice  zone.  One  way  to  work  around  these  large 
rotations  is  to  use  a  method  that  ignores  rotations  in  a 
scene  altogether.  Invariant  moments  is  one  such  method, 
whereby  invariant  moments  reduces  a  scene  to  a  set  of 
numbers,  indices  that  do  not  change  with  rotation. 

The  problem  of  matching  rotated  sea  ice  features  is 
more  complex,  however,  than  the  problem  of  coping  with 
rotational  transformations  because  sea  ice  features  may 
change  with  translation  or  rotation  as  ice  floes  move 
across  open  water  or  ice  fields.  Floes  may  collide  with 
each  other,  thereby  creating  new  features.  Rotations  of 
floes  will  affect  how  features  are  manifested  in  a  SAR 
image  because  SAR  is  sensitive  to  the  aspect  angle  from 
which  a  scene  is  viewed  [2,3].  Translations  of  these 
rotated  floes  might  occur  across  a  changing  background 
of  new  ice,  wind-roughened  ocean,  or  pack  ice.  A  sea  ice 
tracking  algorithm  is  put  to  task  under  these  conditions. 

IPAC  methods  will  handle  a  significant  number  of 
these  rotated  and  translated  ice  features,  as  discussed  in 
the  previous  section.  However,  not  all  ice  features  are 
amenable  to  this  type  of  tracking.  For  example,  ice  floes 
in  the  marginal  ice  zone  commonly  rotate  many  times 
while  drifting  across  open  water.  It  is  for  these  types  of 
ice  features  that  we  consider  invariant  moments.  We  think 
that  invariant  moments  will  also  be  helpful  in  sorting  out 
ambiguities  which  arise  during  application  of  IPAC 
algorithms  to  regions  with  ice  rotation  greater  than  the 
limits  suggested  in  section  2.  . 

Although  a  variety  of  rotation-translation  invariant 
measures  exist,  we  use  those  given  by  Hu  (9).  These 
seven  indices  demonstrate  invariance  under  translation 
and  rotation.  Furthermore,  they  also  demonstrate 
invariance  under  scaling  and  optical  inversion.  Such 
properties— invariance  under  translation,  rotation,  scaling, 
and  mirroring — make  these  indices  well  suited  for  scene 
matching  problems  [1 0].  To  show  why  Hu’s  invariant 
moments  have  these  properties,  we  will  highlight  their 
derivation. 

Hu's  invariant  moments  start  with  the  theorem  stated 
by  PapouBs  in  1968.  .  In  essence,  this  theorem  infers 

that  an  image  represented  by  a  image  intensity  function 


f(x,y)  bounded  in  extent  can  be  reduced  to  a  unique 
sequence  of  double  moments.  This  means  that  an  image 
many  thousands  of  pixels  large  will  compact  down  to  a 
sequence  of  a  few  numbers.  Note  that  low-order  double 
moments  of  images  are  analogous  to  moments  for  objects 
(mass,  center  of  mass,  moments  of  inertia). 

Double  moments  by  themselves  show  no  invariant 
properties.  However,  Hu  derived  invariant  moments  from 
double  moments  with  a  series  of  transformations.  From 
double  moments,  we  obtain  central  moments.  Under  the 
translations  of  coordinates,  central  moments  remain 
invariant.  From  central  moments,  we  obtain  normalized 
central  moments.  Under  scaling,  in  addition  to  translation, 
normalized  central  moments  remain  invariant.  From 
normalized  central  moments  and  with  the  theory  of 
algebraic  invariants,  we  obtain  Hu's  indices.  Under 
proper  and  improper  orthogonal  rotation  and  mirroring,  in 
addition  to  translation  and  scaling,  Hu's  indices  remain 
invariant.  Finally,  since  the  magnitude  of  these  moments 
may  span  many  orders,  the  typical  convention  is  to  take 
the  absolute  value  of  the  log  of  Hu's  indices. 


3.2  Case  Study:  SEASAT  Orbits  1438  and  1481 

While  in  theory,  Hu's  invariant  moments  look 
attractive  as  an  alternative  matching  algorithm  to  IPAC 
methods  in  tracking  some  of  the  more  profoundly  rotated 
sea  ice  features,  in  practice  the  invariant  moments  work 
less  well.  A  significant  problem  with  Hu's  invariant 
moments  has  to  do  with  interference  from  other  floes. 

To  illustrate  how  Hu’s  invariant  moments  work  with 
SAR  images  of  sea  ice,  we  take  several  examples  from 
the  SEASAT  image  pair  from  orbits  1438  and  1481. 

These  particular  examples  come  from  a  marginal  ice  zone 
in  the  Beaufort  Sea. 


Fig.  5.  Example  of  a  rotated  ice  floe,  from  the  marginal  ice  zone  in 
the  Beaufort  Sea.  Image  resolution  was  averaged  to  250  meters.  The 
window  size  is  100x100  pixels,  a)  Reference  floe:  orbit  1438,  October 
5,1978.  b)  Relerence  Floe:  orbit  1481,  October  8, 1978. 

Figure  5a  depicts  an  isolated  floe  (orbit  1438)  which 
will  Serve  as  our  reference  floe.  Figure  5b  depicts  that 
same  floe  which  has  since  drifted  and  rotated  (orbit  1481). 

Rotation  of  this  floe  apparently  did  affect  how  its  features 
are  manifest  in  a  SAR  image:  the  lower  end  of  the  floe 
brightens  while  the  rest  of  the'floe  dims  with  the  rotation. 
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Furthermore,  the  translation  of  this  floe  did  occur  across  a 
changing  background  of  ice  and  wind-roughened  ocean. 
The  ice  floe  in  1481  appears  in  what  looks  like  a  - 
thickening  background  of  ice,  thereby  altering  some  of  the 
floe’s  boundaries:  the  bright  lower  right  hand  corner  of 
Figure  5b  is  interpreted  as  wind-roughened  ocean. 


Fig.  6  Histogram  of  the  absolute  value  of  the  log  of  the  invariant 
moments  corresponding  to  the  reference  floe. 


Figure  7  depicts  two  other  floes  for  comparison.  Both 
floes  were  taken  from  1481.  Figures  8  and  9  illustrate 
how  the  invariant  moments  from  these  floes  compare  with 
the  invariant  moment^  from  the  reference  floe  in  1438. 


Fig.  8.  Histogram  of  the  absolute  value  of  the  log  of  the  invariant 
moments  corresponding  to  Floe  A.  The  invariant  moments  from  the 
reference  floe  1438  are  included  for  comparison.  Invariant  moments 
correctly  distinguish  between  this  and  the  reference  floe. 


Figure  6  shows  a  histogram  of  the  invariant  moments 
corresponding  to  Figure  5.  This  histogram  illustrates  two 
points.  The  first  point  is  that  the  invariant  moments  span 
about  seventeen  orders  of  magnitude,  ranging  from  10’^ 
to  KT1^  This  large  range  of  values  is  common  with  Hu’s 
invariant  moments.  In  general,  the  first  invariant  moment 
has  the  lowest  order  of  magnitude,  with  the  order  of 
magnitude  increasing  with  order  of  invariant  moment. 

The  second  point  is  that  the  invariant  moments  match  well 
at  the  first,  third,  fourth,  and  seventh  moments.  As  a 
rough  indicator,  these  matches  between  invariant 
moments  do  point  to  a  correct  match  between  scenes. 


Figure  8  illustrates  how  dissimilar  scenes 
might  appear  in  terms  of  invariant  moments.  With  the 
exception  of  the  first  invariant  moment,  all  the  remaining 
invariant  moments  of  floe  A  fell  below  those  of  our 
reference  floe.  Other  dissimilar  scenes  would  have  a 
invariant  moment  signature  different  from  the  reference 
signature. 


Invariant  Moment 


Fig.  7.  Others  floes  from  orbit  1481 .  Image  resolution  was 
averaged  to  250  meters.  The  window  size  is  100x100  pixels,  a)  Floe 
A.  Invariant  moments  correctly  distinguish  between  this  and  the 
reference  floe  1438.  b)  Floe  B.  Invariant  moments  yield  ambiguous 
results  between  this  and  the  reference  floe  1438. 


Fig.  9.  Histogram  of  the  absolute  value  of  the  log  of  the  invariant 
moments  corresponding  to  Floe  B.  The  invariant  moments  from  the 
reference  floe  1438  are  included  for  comparison.  Invariant  moments 
yield  ambiguous  results  between  this  and  the  reference  floe. 
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Figure  9  illustrates  how  ambiguous  scenes  might 
appear.  Even  though  the  two  scenes  look  dissimilar,  the 
invariant  moments  of  floe  B  suggest  that  the  reference  floe 
and  floe  B  are  matched.  This  ambiguity  could  occur 
when  the  background,  in  this  case  a  part  of  another  floe 
above  floe  B,  interferes  with  determining  a  match.  This 
problem  would  occur  with  correlation  methods  as  well. 
Other  ambiguities  would  also  occur  should  a  smaller  flow 
look  similar  to  the  reference  floe:  Hu’s  invariant  moments 
are  also  scale  invariant. 

This  work  in  invariant  moments  illustrates  what  kind 
of  invariant  moment  signatures  we  might  get  from 
applying  Hu's  invariant  moments.  It  also  indicates 
possible  sources  of  ambigous  matches.  For  further  work, 
we  need  to  investigate  methods  that  reduce  the  likelihood 
of  ambiguous  matches.  For  example,  some 
improvements  to  invariant  moment  methods  for  scene 
matching  problems  are  given  by  Goshtasby  (11].  Such 
methods  need  to  be  applied  to  Hu's  invariant  moments 
before  they  can  be  considered  as  replacements  for  block 
correlation  measures. 


4.  FEATURE  TRACKING  WHEN  SEA  ICE 
ROTATION  IS  PRESENT 

4.1.  Introduction  to  the  Method 

The  feature  tracking  algorithm  approaches  the 
problem  of  sea  ice  tracking  by  relying  on  selected  sea  ice 
features  rather  than  the  entire  content  of  the  image,  as  in 
the  IPAC  method  of  section  2.  There  are  a  wide  variety  of 
features  in  sea  ice  which  one  could  track,  including 
pressure  ridges,  isolated  floes,  and  floe-lead  boundaries. 
Initially  we  have  chosen  floe-lead  boundaries  as  the 
physical  feature  we  will  track  because  they  are  usually  the 
most  prominent,  distinct,  and  high  contrast  featurfes  in 
SAR  images  of  sea  ice.  The  movement  of  these  features 
from  one  image  to  the  next  provides  the  information  from 
which  the  sea  ice  velocity  field  is  estimated.  Feature 
patterns  must  be  very  distinctive,  if  not  unique,  objects. 

In  the  feature  tracking  method  each  image  is 
transformed  into  a  set  of  features  described  by  parameters 
independent  of  the  feature  location  and  orientation.  The 
collection  of  several  hundred  features,  each  described  by 
several  hundred  numbers,  reduces  the  information 
needed  to  characterize  a  100  x  100  km  image  by  roughly 
two  orders  of  magnitude.  The  search  space  is 
correspondingly  reduced.  Once  a  feature  pattern  is 
recognized  in  both  members  of  a  pair  of  images,  sea  ice 
displacement  and  hence  velocity  information  can  be 
derived. 

In  the  following  paragraph  we  outline  our  feature 
tracking  algorithm;  more  detailed  information  is  given  in 
(5,6].  Given  a  pair  of  SAR  images  of  the  same  geographic 
area,  each  image  is  first  averaged  to  a  resolution  of  about 
250  m.  Ice  and  lead  are  then  classified  by  a  simple 
intensity  threshold  and  floe-lead  boundaries  identified. 


Pixels  along  floe-lead  boundaries  are  classified  and 
assembled  into  segments.  These  boundary  segments  are 
characterized  and  thus  form  the  features  which  are 
tracked  from  one  SAR  image  to  another  collected  three 
days  later.  Several  tracking  passes  are  made  through  the 
collections  of  features  in  each  image.  The  high 
confidence  ice  displacement  estimates,  corresponding  to 
high  correlation  between  the  characteristics  of  matching 
features  in  the  two  images,  are  used  to  guide  the  search 
for  other  matching  features  and  thus  more  ice 
displacement  estimates. 


4.2  Capabilities  of  Feature  Tracking  with  Regard 

to  Rotation 

The  feature  tracking  algorithm  is  inherently  robust 
with  respect  to  rotation  because  the  features  are 
characterized  with  respect  to  a  coordinate  system  tied  to 
the  feature  itself  and  not  with  respect  to  geographic  or 
image  coordinates.  To  understand  this  we  consider  the 
characterization  process  in  more  detail.  After  pixels  along 
floe-lead  boundaries  are  classified  and  assembled  into 
segments  these  boundary  segment  features  are 
characterized  by  their  deviation  from  a  straight  line.  First, 
a  'least  squares'  criterion  straight  line  is  fit  to  a  given 
boundary  segment.  The  perpendicular  distance  from  this 
straight  line  to  the  segment  is  calculated  at  regular 
intervals  along  the  straight  line.  This  set  of  deviations 
characterizes  the  shape  of  a  boundary  segment  as 
illustrated  in  Figure  10.  This  method  of  characterizing  a 
boundary  segment  feature  is  invariant  under  rotation  of 
the  feature.  There  are  limitations,  however,  as  noted  in 
section  4.3. 

Recognition  of  matching  features  is  done  simply  by 
comparing  the  sets  of  deviation  measurements  for  two 
candidate  features.  A  variance  normalized, 
cross-correlation  coefficient  (C)  is  used  to  evaluate  how 
closely  two  candidate  features  match.  Often  segments 
which  we  want  to  compare  are  of  different  lengths 
because  of  partial  breakup  of  a  boundary  segment  over 
the  time  interval  between  images.  So  in  our  search  we 
include  features  of  different  lengths.  During  the 
comparison  process  we  slide  the  shorter  feature  along  the 
longer  feature  with  C  calculated  for  all  placements  of  the 
smaller  feature  with  respect  to  the  larger.  This  allows  for 
larger  segments  to  break  into  smaller  ones  and  still  be 
found.  A  successful  match  reveals  ice  displacement. 

Because  this  method  of  characterizing  a  boundary 
segment  feature  is  invariant  under  rotation  of  the  feature 
one  expects  that,  in  principle,  a  feature  which  has 
undergone  rotation  will  be  as  easily  found  as  one  which 
has  not.  In  the  one  pair  of  SAR  sea  ice  images  so  far 
examined  (5,6]  features  were  successfully  tracked  using 
this  algorithm.  However,  the  SAR  scenes  contained 
central  pack  ice  and  the  average  rotation  of  the  sea  ice 
was  only  about  5°. 
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Fig.  10.  Boundary  segment  lilted  by  a  straight  line  and 
characterized  by  deviation  from  the  lilted  line 


4.3  Limitations  of  Feature  Tracking  with  Regard 
to  Rotation 


left  and  right  edges  will  be  different  because  a  pixel  which 
was  type  5  in  the  first  image  might  become  type  4  or  6  m 
the  second  image  of  the  pair.  Under  these  circumstances 
a  rotated  feature  becomes  more  difficult  to  track.  If  the 
majority  of  sea  ice  in  a  region  rotates  approximately  the 
same  amount,  then  the  algorithm  will  adapt  to  the  new, 
dominant,  lead  direction.  Tracking  in  the  presence  of 
rotation  should  then  proceed  as  in  the  case  of  no  rotation. 


Although  our  feature  tracking  algorithm  has 
advantages  in  terms  of  following  rotating  ice,  there  are 
also  limitations  so  far  as  the  current  algorithm  [5,6]  is 
concerned.  The  limitation  arises,  not  in  the 
characterization  of  a  feature,  but  in  choosing  what  pixels 
constitute  the  feature.  Features  are  constructed  by  first 
classifying  floe-lead  boundary  pixels  and  then  assembling 
them  into  contiguous  sequences  of  pixels  having  only 
certain  allowed  classifications.  The  classification  system 
uses  surrounding  pixels  for  reference  and  thus  is  not 
completely  independent  of  image  coordinates.  This 
process  allows  two  sides  of  a  lead  to  be  separate  features 
and  thus  move  independently,  i.e.  it  allows  leads  to  open 
and  close.  Thus  the  lead-floe  boundary  surrounding  a 
lead  may  be  chopped  into  features  (for  tracking)  in  two 
different  ways  depending  on  how  it  is  oriented  in  the  first 
and  second  images  of  an  image  pair. 
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Fig.  1 1 .  Diagram  illustrating  edge  pixel  classification  scheme.  The 
filled  circles  represent  water  pixels  and  the  white  circles  Ice  pixels. 

This  problem  can  be  illustrated  by  considering  the 
diagram  of  Figure  1 1 .  Leads  tend  to  be  cracks  in  ice  and 
thus  linear  rather  than  circular.  If  the  majority  of  the  leads 
in  an  image  are  linear,  running  vertically,  then  our 
algorithm  [5,6]  will  assemble  edge  features  as  follows: 

Left  Edges:  edge  pixel  types  2,3,4  and  associated 
types  1  &  5 

Right  Edges:  edge  pixel  types  6,7,8  and  associated 
types  1  &  5.. 

If,  in  the  second  image  of  a  pair,  a  particular  lead  has 
rotated  by  say  45°;  then  some  of  the  pixels  making  up  the 


4.4  Improvements  of  Feature  Tracking  with 

Regard  to  Rotation 

Feature  tracking  is  a  broad  catagory  of  pattern 
recognition  techniques  for  application  to  sea  ice.  As 
discussed  above,  it  has  inherent  advantages  for  ice 
tracking  in  the  presence  of  rotation.  We  also  discussed  a 
limitation  which  arises  from  the  fact  that  leads  often  open 
and  close  and  thus  the  dominant  edges  must  be  tracked 
independently.  Other  sea  ice  features  do  not  have  this 
limitation.  For  example,  pressure  ridge  patterns  are  an 
evident  feature  in  figure  1  above.  Although  pressure 
ridges  are  less  distinct  features  than  lead-floe  boundaries, 
they  are  more  stable  since  they  are  not  affected  by 
opening  and  closing  of  leads.  Thus  one  would  expect  that 
incorporating  pressure  ridge  patterns  into  a  feature 
tracking  algorithm  would  improve  performance, 
particularly  in  the  presence  of  large  and  spatially  variable 
ice  rotation. 


5.  CONCLUSIONS 

We  summarize  our  major  conclusions  as  follows: 

1.  The  image  pyramid  area  correlation  (IPAC) 
method,  also  known  as  hierarchical  correlation,  has 
successfully  performed  automated  sea  ice  tracking;  but  is 
vulnerable  to  uncertainties  if  the  ice  rotates  more  than 
about  10°  to  15°  between  SAR  observations. 

2.  Window  size  is  important  in  making  block 
cross-correlation  comparisons  between  SAR  images.  We 
find  that  because  ice  motion  is  nonuniform  window  size 
choice  is  a  tradeoff  between  better  placement  accuracy 
with  smaller  windows  and  lower  levels  of  spurious 
sidelobe  response  with  larger  windows. 

3.  Rotation  invariant  correlation  methods  can 
successfully  track  rotating  ice  floes  which  are  isolated. 
This  technique  may  also  be  useful  is  resolving 
ambiguities  in  IPAC  displacement  estimates  which  are 
due  to  ice  rotation  (see  item  1  above). 

4.  Feature  tracking  is  inheirantly  robust  with 
respect  to  tracking  rotating  sea  ice,  but  has  some 
limitations  when  floe-lead  boundaries  are  used  as 
features.  The  use  of  pressure  ridge  features  is  indicated. 
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5.  We  conclude  that  an  automated  sea  ice 
tracking  algorithm  will  need  a  variety  of  techniques  to 
operate  successfully  over  a  wide  range  of  sea  ice 
conditions. 
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OBSERVATION  OF  SEA  ICE  MOTION  AND  DEFORMATION 
USING  SYNTHETIC  APERTURE  RADAR: 
AUTOMATED  ANALYSIS  ALGORITHMS 
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Abstract  --  Synthetic  aperture  radar  (SAR)  provides  an  excellent  means  of  observing  the 
movement  and  distortion  of  sea  ice  over  large  temporal  and  spatial  scales.  SAR  observations  are 
not  affected  by  clouds  or  lack  of  sunlight.  The  European  Space  Agency's  ERS-1  satellite  as  well  as 
others  planned  to  follow  will  carry  SARs  over  the  polar  regions  beginning  in  1990.  A  key 
component  in  utilization  of  these  SAR  data  is  an  automated  scheme  for  extracting  the  sea  ice 
velocity  field  from  a  time  sequence  of  SAR  images  of  the  same  geographical  region.  In  this  paper 
we  provide  an  overview  of  the  problem  of  automated  sea  ice  tracking  using  SAR  images.  To  begin 
we  briefly  describe  two  basic  techniques  for  automated  sea  ice  tracking:  image  pyramid  area 
correlation  (hierarchical  correlation)  and  feature  tracking.  Examples  of  the  application  of  these 
techniques  to  SEASAT  SAR  sea  ice  images  are  shown.  The  results  compare  well  with  each  other 
and  with  manually  tracked  estimates  of  the  ice  velocity  field.  We  note  the  problems  presented  to 
automated  tracking  algorithms  by  sea  ice  which  rotates  or  shears,  discussing  the  dimensions  of 
the  problem  as  well  as  solutions  via  the  use  of  rotation  invarient  correlation  and  feature 
tracking.  Using  the  ice  velocity  field  estimated  from  SAR  images  one  can  reconstruct  one  sea  ice 
image  from  the  other  member  of  the  pair.  Comparing  the  reconstructed  image  with  the  observed 
one,  errors  in  the  estimated  velocity  field  can  be  recognized  and  a  useful  probable  error  display 
created  automatically  to  accompany  ice  velocity  estimates  provided  to  the  user.  We  suggest  that 
this  error  display  may  be  useful  in  segmenting  the  sea  ice  observed  into  regions  which  move  as 
rigid  plates  and  regions  of  significant  ice  velocity  shear  and  distortion.  We  then  consider  how  the 
several  methods  discussed  can  fit  together  to  form  a  robust  algorithm  for  application  to  satellite  (ERS-1) 
SAR  sea  ice  data  to  be  collected  at  the  Alaskan  SAR  FAcility  (ASF)  at  the  University  of  Alaska,  Fairbanks. 
Finally,  we  discuss  the  importance  of  displaying  the  ice  velocity  field  properly  and  summarize  conclusions. 


I.  INTRODUCTION 

The  importance  of  sea  ice  remote  sensing  arises  from  the  high  spatial  and 
temporal  variability  of  sea  ice,  its  wide  geographic  extent  and  its  impact  on  climate, 
arctic  off-shore  engineering,  transportation  and  military  operations.  Sea  ice 
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dramatically  changes  the  air-sea  energy  balance  and  energy  transport.  Since  sea  ice 
can  cover  significant  areas  of  the  Arctic  and  Southern  Oceans,  it  has  an  important 
impact  on  both  polar  and  global  climate. 

The  high  variability  of  sea  ice  coupled  with  its  large  extent  makes  remote  sensing 
techniques  a  necessity  for  large-scale  studies.  Since  polar  regions  are  often  cloud 
covered  and  without  sunlight  for  long  periods,  the  ability  of  synthetic  aperture  radar 
(SAR)  to  provide  areally  extensive,  high  resolution  images  through  clouds  and  without 
sunlight  gives  it  a  central  role  in  sea  ice  remote  sensing.  SAR  provides  good 
information  on  sea  ice  extent,  movement  and  deformation,  internal  geometry  (structure 
of  floes  and  leads)  and  surface  roughness  (Weller,  et  a!.,  1983).  Some  useful 
information  on  sea  ice  types  and  their  physical  properties  can  also  be  obtained  from 
SAR  images. 

The  process  of  retrieving  motion  information  from  a  sequence  of  SAR  images  (of 
the  same  geographic  area)  revolves  around  two  principal  items,  namely: 

1.  Identifying  regions  in  a  sea  ice  field  to  use  as  'tie  points'  for  ice 

floe  tracking. 

2.  Finding  the  same  tie  points  in  subsequent  images  of  the  sequence. 

Once  the  ice  displacement  is  known  the  ice  velocity  can  easily  be  calculated  since  the 
elapsed  time  between  the  SAR  images  is  known. 

This  process  has  been  carried  out  manually  on  SEAS  AT  SAR  data  by  Thorndike 
and  Rothrock  [see  Weller  et  al.,  1983]  using  a  split-screen  video  display  and  an 
experienced  operator.  Ice  velocity  is  computed  automatically  once  the  locations  of  tie 
points  are  marked  in  a  pair  of  images  on  the  screen.  Further  work  on  satellite  imagry 
using  manual  techniques  has  been  done  by  Carsey  et  al.  (1983)  and  by  Curlander, 

Holt  and  Hussey  (1987).  However,  this  process  is  too  slow  and  expensive  for  real  time 
and/or  extensive  analysis  of  images  provided  by  long-term  satellite  SARs.  Two  and 
possibly  three  such  satellite  SARs  are  scheduled  to  be  in  orbit  by  the  1990s,  namely  the 
European  Space  Agency  ERS-1 ,  the  Japanese  ERS-1  and  the  joint  Canadian-US-UK 
RADARSAT.  Real  aperture  radar  images  at  low  resolution  are  apparently  available 
from  the  Soviet  Kosmos  1500  series  satellites.  Hence,  the  development  of  automated 
techniques  for  sea  ice  tracking  is  a  primary  requirement  for  successful  use  of  satellite 
SAR  images  in  sea  ice  research. 

In  this  paper  we  present  an  overview  of  the  ice  tracking  problem,  discussing  the 
several  component  methods  used  in  automated  sea  ice  tracking,  detection  of  errors,  the 
structure  of  composite  algorithms  for  ice  tracking  under  a  wide  range  of  sea  ice 
conditions  and  finally  the  display  of  sea  ice  velocity  field  estimates.  Sections  II,  III  and 
IV  cover  image  pyramid  area  correlation  and  feature  tracking  methods  as  well  as 
methods  for  coping  with  situations  in  which  the  sea  ice  rotates  between  observations. 
Section  V  discusses  the  important  problem  of  error  detection  using  an  image 
reconstruction  method.  The  synthesis  of  a  composite  algorithm  is  covered  in  section  VI. 
Section  VII  deals  with  the  display  of  velocity  field  results  so  as  to  aid  the  user  with 
interpretation.  Conclusions  are  presented  in  section  VIII. 
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II.  IMAGE  PYRAMID  AREA  CORRELATION  METHOD 


A.  Description  of  Method 

This  method,  also  known  as  'hierarchical  correlation'  or  'nested  correlation', 
comes  from  the  field  of  pattern  recognition  and  was  first  applied  to  the  tracking  of  sea 
ice  by  Fily  and  Rothrock  (1986).  Cross-correlation  of  small  blocks  of  pixels  from  a  pair 
of  cospatial  images  is  used  to  find  matching  scene  points.  This  is  made  computationally 
efficient  by  first  building  from  each  primary,  high  resolution  image  a  sequence  (pyramid) 
of  derived  images  of  progressively  coarser  resolution,  but  correspondingly  fewer  data 
points.  Such  a  pair  of  image  pyramids  is  shown  in  Fig.  1 .  Block  cross  correlation  in  the 
coarsest  resolution  images  of  the  pyramid  is  computationally  fast  and  builds  a  crude 
model  of  the  displacement  field  of  the  sea  ice.  This  crude  model  is  then  used  to  guide 
searches  in  the  pyramid  of  progressively  higher  resolution  images.  As  one  moves 
down  the  pyramid  to  the  highest  resolution  image,  searches  are  guided  from  previous 
results  and  the  displacement  field  is  continually  refined. 


Fig.  1 .  Image  pyramids  for  two  SEASAT  SAR  images  of  the  same  central  pack  ice  area  in  the  Beaufort 
Sea.  The  images  were  collected  three  days  apart  on  Oct.  5  (top)  &  8  (bottom),  1978  on  Seasat  orbits  1439 
and  1482,  respectively. 


An  image  pyramid,  as  shown  in  Fig.  1 ,  can  be  viewed  as  a  set  of  images  in  a  stack. 
At  the  bottom  of  the  stack  is  the  original  image  and  at  each  of  the  following  levels  of  the 
stack  is  a  lower  resolution  image  derived  from  the  level  below  it  by  first  applying  a  low 
pass  filter  and  then  resampling  the  image.  The  result  is  a  pyramid  in  which,  as  one 
traverses  from  bottom  to  top,  the  number  of  pixels  in  the  images  decreases 
geometrically. 

The  following  steps  are  needed  for  effective  correlation  in  a  pyramid  structure.  To 
begin,  pyramids  (1 ,  first  in  time  and  2,  second  in  time)  are  constructed  as  in  Fig.  1 ;  one 
for  each  of  the  two  images.  Then  the  following  algorithm  is  used: 

1 .  For  each  I  x  I  pixel  sample  block  in  the  coarsest  level  of  pyramid  1 ,  calculate  a  displacement 

vector  by  finding  the  global  maximum  of  the  correlation  with  the  coarsest  level  of  Pyramid  2. 

Generate  a  confidence  measure  for  the  displacements. 

2.  For  each  finer  resolution  level  of  the  pyramid  do  the  following: 

a.  Smooth  the  displacement  vector  field  for  the  previous  (coarser)  pyramid  level  in  order  to 
remove  gross  errors.  This  helps  prevent  errors  in  the  top  levels  of  the  pyramid  propagating 
to  the  lower  levels. 

b.  Interpolate  the  displacement  vector  field  calculated  at  the  previous  pyramid  level  to  provide 
displacement  estimates  for  each  I  x  I  pixel  sample  block  of  the  current  level  of  pyramid  1. 

c.  For  each  I  x  I  pixel  sample  block  of  the  current  level  of  Pyramid  1  find  the  local  maximum  of 
the  correlation  with  the  current  level  of  Pyramid  2  by  searching  in  an  M  x  M  window 
centered  at  the  displacement  estimate  of  the  current  level  of  Pyramid  2. 

The  algorithm  starts  at  the  top,  coarsest  pyramid  level  and  proceeds  'top-down'  one 
level  at  a  time.  At  the  coarsest  level,  the  area  to  be  searched  in  the  second  image  is  not 
constrained,  i.e.  M  is  the  dimension  of  the  entire  image.  It  is  very  important  to  have  valid 
estimates  at  this  stage  since  the  following  stages  are  dependent  on  the  first.  At  each 
successive  level  in  the  pyramid,  the  estimated  location  of  the  largest  correlation  found  at 
the  previous,  coarser  pyramid  level  constrains  the  locations  searched.  These  location 
constraints  greatly  decrease  the  total  number  of  locations  that  must  be  searched.  The 
overall  gain  in  speed  over  brute-force  correlation  is  between  100  and  1000. 

The  selection  of  the  size  (I)  of  the  correlation  sample  block  involves  a  trade-off 
between  accurate  determination  of  the  location  of  correlation  maximum  (small  I)  and 
avoidance  of  false  correlation  maxima  (large  I).  Vesecky  et  al.  (1987)  describe  a  cross 
correlation  technique  for  choosing  I.  The  choice  of  I  is  dependent  on  image  content. 

For  the  sea  ice  image  of  Fig.  1 ,  an  I  corresponding  to  about  a  4x4  km  block  size  was 
judged  to  be  near  optimum,  i.e.,  I  =  1 6  for  250  m  size  pixels.  This  size  is  of  the  order  of 
the  size  of  the  average  ice  floe. 

The  quantity  M  is  selected  to  correspond  to  a  search  area  of  the  order  of  the 
physical  size  of  the  correlation  sample  block  in  the  coarser  pyramid  level  just 
processed.  We  have  taken  M  such  that  the  search  area  is  4  times  larger  than  the 
previous  (coarser)  level  sample  block  size.  Further  information  on  the  IPAC  method  is 
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given  by  Fily  and  Rothrock  (1986),  Vesecky,  et  a!.  (1986, 1987  &  1988)  and  Samadani 
(1987). 

B.  Confidence  Estimates 

There  are  two  parts  of  the  pyramid  algorithm  where  confidence  measures  are 
useful.  For  the  final  vector  output,  confidence  measures  can  be  used  to  prune  or  to  at 
least  highlight  suspect  vectors.  In  addition,  during  the  interpolations  from  a  coarser 
pyramid  level  to  a  finer  pyramid  level,  the  confidence  measures  can  be  used  as  a 
weighting  factor  for  the  interpolation  (Mostafavi  and  Smith, 1978).  The  first  technique 
was  used  in  obtaining  the  results  of  Fig.  2.  It  is  shown  by  Vesecky  et  al.  (1988)  that 
there  is  a  correspondence  between  areas  of  low  confidence  and  errors  in  the 
displacement  field  estimated  by  the  IPAC  method. 

C.  Smoothing  Techniques 

In  the  "real"  world  of  sea  ice  remote  sensing  using  SAR  images  one 
encounters  noise,  grey  level  variation  across  a  given  image  and  between  images  as 
well  as  gross  geometric  distortions.  These  problems  result  in  errors  when  block 
correlation  is  used  for  finding  ice  displacement  vectors.  In  view  of  these  problems,  the 
assumption  is  usually  made  that  the  displacement  field  is  smooth  for  most  areas  of  the 
image  (Borodachev.l  984).  Smoothing  techniques  may  be  used  to  remove  the  gross 
errors  in  the  displacement  field.  These  techniques  are  especially  important  for  pyramid 
correlation  because  errors  made  at  the  coarser  resolutions  propagate  to  the  finer  levels 
due  to  the  constraints  on  search  locations  at  the  finer  levels.  One  must  be  careful  when 
using  smoothing  techniques,  however,  in  areas  where  there  is  high  shear. 

In  the  algorithm  described  in  section  II.A  above,  smoothing  of  the  estimated  ice 
velocity  field  was  done  at  each  pyramid  level  by  the  following  procedure.  The  coarse 
level  vectors  were  replaced  by  the  median  of  the  surrounding  high  confidence  vectors. 
Bilinear  interpolation  was  then  used  to  fill  in  vectors  for  the  next  higher  resolution  level 
in  the  pyramid.  At  the  highest  resolution  level  of  the  pyramid  this  procedure  produced  a 
preliminary  result  leading  to  Fig.  2. 

To  achieve  the  result  of  Fig.  2  we  applied  a  different  smoothing  technique  to  the 
preliminary  result  above  by  extending  the  controlled-continuity  splines  method  of 
Terzopoulos  (1985)  for  fitting  vector  fields.  This  more  versatile  interpolation  technique 
makes  systematic  use  of  the  confidence  measure  and  minimizes  the  weighted  square 
deviation  from  the  displacement  data  while  at  the  same  time  minimizing  the  integral  of 
the  second  spatial  derivatives.  Samadani  (1 986)  gives  further  discussion  of  this 
smoothing  technique. 

D.  Results  for  SEASAT  SAR  Image  Pair 

The  pyramid  algorithm  described  above  was  applied  to  two  averaged  (250 
meter  pixels)  SEASAT  images  of  central  arctic  pack  sea  ice  -  Rev.  1439  (October  5, 
1978)  and  Rev.  1482  (October  8,  1978),  Fig.  1,  by  Vesecky  et  al.  (1986  &  1988).  Fig.  2 
compares  these  displacement  vectors  with  a  computer  assisted  manual  analysis  by 
Thorndike  and  Rothrock  [see  Weller  et  al.,  1983];  it  is  seen  that  the  method  performs 
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very  well.  However,  there  are  areas  where  the  automated  algorithm  has  generated 
erroneous  velocity  vectors.  The  problem  areas  at  the  edges  of  the  figure  correspond  to 
ice  features  flowing  into  or  out  of  the  finite  area  of  these  images.  Further  discussion  of 
errors  is  included  in  section  V.  The  results  of  Fig.  2  can  be  used  to  calculate  the 
divergence  and  curl  of  the  velocity  vector  field  as  done  by  Vesecky  et  al.  (1988).  These 
fields  in  turn  can  be  used  to  investigate  sources,  sinks  and  vorticity  in  the  ice  flow  as 
well  as  appearance  and  disappearance  of  open  water. 


Fig  2.  (left)  Sea  ice  velocity  vectors  generated  by  the  image  pyramid  area  correlation  technique  (IPAC) 
applied  to  the  SAR  images  of  Fig.  1.  (right)  sea  ice  velocity  vectors  generated  from  the  same  image  data  of 
Fig.  1  by  Thorndike  and  Rothrock  using  a  computer  assisted  manual  technique.  The  IPAC  velocity 
estimates  above  have  been  improved  by  pruning  of  low  confidence  estimates  and  using  the 
controlled-continuity  spline  technique  for  interpolation. 


III.  FEATURE  TRACKING 
A.  Introduction  to  the  Method 

Feature  tracking  approaches  the  problem  of  sea  ice  motion  estimates 
somewhat  differently  by  relying  on  selected  sea  ice  features  rather  than  the  content  of 
the  entire  image.  There  are  a  wide  variety  of  features  in  sea  ice  which  one  could  track, 
including  pressure  ridges,  isolated  floes  and  floe-lead  boundaries.  The  results  of  Fig.  3 
below  rely  on  floe-lead  boundaries  as  the  physical  feature  to  be  tracked  because  they 
are  usually  the  most  prominent,  distinct  and  high  contrast  features  in  SAR  images  of  sea 
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ice.  The  movement  of  these  features  from  one  image  to  the  next  provides  the  'tie  points' 
from  which  the  sea  ice  velocity  field  is  estimated.  Feature  patterns  must  be  very 
distinctive,  if  not  unique,  objects.  They  are  characterized  by  information  which  is 
independent  of  the  object  location  and  orientation,  though  this  latter  information  is,  of 
course,  attached  to  them.  The  collection  of  several  hundred  features,  each  described 
by  several  hundred  numbers,  reduces  the  information  needed  to  characterize  a  100  x 
100  km  SAR  image  by  roughly  two  orders  of  magnitude.  The  search  space  is 
correspondingly  reduced. 

Once  a  feature  pattern  is  recognized  in  both  members  of  a  pair  of  images,  sea  ice 
displacement  and  hence  velocity  information  can  be  derived.  Since  we  know  the 
movement  of  every  pixel  in  the  feature,  we  can  derive  not  just  one,  but  many  velocity 
vectors.  Further,  we  know  the  rotation  as  well  as  the  displacement  of  a  feature,  since 
we  keep  track  of  a  feature's  geographic  orientation  in  both  images  of  the  pair.  The 
results  shown  in  Fig.  3  below  are  for  floe-lead  boundaries;  however  the  method  can 
also  be  applied  to  other  features  such  as  pressure  ridges  or  isolated  floes. 

B.  Description  of  a  Feature  Tracking  Algorithm 

Here  we  briefly  describe  the  feature  tracking  algorithm  of  Vesecky  et  al.  (1988) 
which  uses  floe-lead  boundaries  as  the  features  to  be  tracked.  Each  of  the  pair  of  SAR 
images  is  first  averaged  to  a  resolution  of  about  250  m.  Ice  and  lead  are  then  classified 
by  a  simple  intensity  threshold  and  floe-lead  boundaries  identified.  Pixels  along 
floe-lead  boundaries  are  classified  and  assembled  into  segments.  These  boundary 
segments  are  characterized  and  thus  form  the  features  which  are  tracked  from  one  SAR 
image  to  another  collected  three  days  later.  Several  tracking  passes  are  made  through 
the  collections  of  features  in  each  image.  The  high  confidence  ice  displacement 
estimates,  corresponding  to  high  correlation  between  the  characteristics  of  matching 
features  in  the  two  images,  are  used  to  guide  the  search  for  other  matching  features  and 
thus  more  ice  displacement  estimates. 

First,  both  SAR  images  of  a  pair  are  averaged  to  increase  the  pixel  size  from  the 
original  25x25  m  size  to  250x250  m  size.  Averaging  reduces  the  speckle  variance 
inherent  in  SAR  images  and  make  computation  more  tractable  since  there  are  a  factor 
of  1 00  fewer  pixels  to  deal  with  after  averaging.  Since  SAR  pixels  cannot  be  located 
geographically  to  better  than  about  200  m,  the  averaging  process  results  in  little  loss  of 
velocity  measurement  accuracy. 

A  histogram  of  pixel  intensity  is  formed  and  a  simple  pixel  intensity  threshold 
chosen  to  distinguish  between  ice  floes  and  leads.  Fily  and  Rothrock  (1986)  have 
shown  that  this  technique  is  viable  for  central  pack  ice.  Thus  the  SAR  image  pair  can 
be  transformed  into  a  pair  of  images  with  binary  pixels,  i.e.  1  for  ice  floe  pixels  and  0  for 
lead  water  pixels  (including  new  ice).  Floe-lead  boundary  pixels  are  identified  simply 
by  finding  water  pixels  adjacent  to  ice  pixels. 

Next,  the  ice-water  boundary  pixels  identified  above  are  classified  into  eight  types 
according  to  their  surroundings  to  determine  the  orientation  of  the  floe-lead  boundary 
edge.  The  four  nearest  neighbor  pixels  are  considered  to  find  if  there  is  water  above, 
below,  to  right  or  to  left  of  a  given  boundary  pixel.  The  classified  pixels  are  then 
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assembled  into  boundary  segments  which  are  the  features  to  be  tracked.  The  assembly 
scheme  begins  by  finding  the  dominant  boundary  direction.  In  Fig.  1  the  dominant 
direction  is  vertical,  so  the  boundary  segments  are  divided  into  two  groups, 
corresponding  to  left  and  right  edges  of  the  leads  in  Fig.  1. 

The  goal,  then,  would  be  to  match  left  edges  in  one  image  with  left  edges  in  the 
other  image  and,  independently,  right  edges  in  the  two  images.  This  allows  the  ice 
floes  along  the  left  side  of  a  lead  to  slip  with  respect  to  ice  floes  along  the  right  side. 

Thus  shear  motion,  which  we  would  expect  across  leads  in  central  pack  ice,  is  built  into 
the  algorithm  for  this  type  of  sea  ice. 

Boundary  segment  features  are  characterized  by  their  deviation  from  a  straight  line. 
First,  a  'least  squares'  criterion  straight  line  is  fit  to  a  given  boundary  segment 
assembled  as  above.  The  distance  from  this  straight  line  to  the  segment  is  calculated  at 
regular  intervals  along  the  straight  line.  This  set  of  deviations  characterizes  the  shape 
of  a  boundary  segment.  Note  that  this  method  of  characterizing  a  boundary  segment 
feature  is  invariant  under  rotation  of  the  feature. 

Recognition  of  matching  features  is  done  simply  by  comparing  the  sets  of  deviation 
measurements  for  two  candidate  features.  A  variance  normalized,  cross-correlation 
coefficient  (C)  was  used  to  evaluate  how  closely  two  candidate  features  matched.  Often 
segments  which  we  want  to  compare  are  of  different  lengths  because  of  partial  breakup 
of  a  boundary  segment  over  the  three  days  between  images.  So  in  the  search,  features 
of  different  lengths  were  included.  During  the  comparison  process  the  shorter  features 
slide  along  the  longer  features  with  C  calculated  for  all  placements  of  the  smaller 
feature  with  respect  to  the  larger.  This  allows  for  larger  segments  to  break  into  smaller 
ones  and  still  be  found. 

The  search  and  matching  procedure  is  an  iterative  one  and  has  features  which 
remind  one  of  the  image  pyramid  method  of  section  II  above.  During  the  several 
iterations  of  the  search  we  use  the  high  confidence  ice  displacement  estimates  of  the 
initial  iterations  to  guide  the  searches  of  subsequent  iterations.  In  the  first  iteration  only 
long  boundary  segments  are  used  since  the  larger  number  of  pixels  gives  a  more  stable 
estimate  of  C  and  a  correct  match  yields  ice  displacement  estimates  over  a  larger 
region.  The  search  is  limited  to  a  displacement  of  60  km--20  km  per  day  for  the  three 
days  between  the  images.  This  search  area  could  be  significantly  reduced  if 
environmental  data  such  as  local  winds  were  available  to  indicate  probable  ice 
movement. 


C.  Results  for  SEASAT  SAR  Image  Pair  in  the  Beaufort  Sea 

Fig.  3  shows  the  results  of  applying  the  boundary  segment  feature  tracking 
algorithm  to  the  pair  of  sea  ice  images  shown  in  Fig.  1 .  Only  the  highest  resolution, 
250x250  m  pixel  size,  images  are  used.  Boundary  segments  constructed  from  these 
images  vary  in  length  from  4  to  163  pixels  or  1  to  41  km  in  length  and  have  a  dominant 
direction  which  is  vertical,  i.e.  North-South,  in  Fig.  1.  The  resulting  sets  of  boundary 
segment  features  in  the  reference  (first  in  time)  and  test  (second  in  time)  images  are 
shown  in  left  and  middle  panels  of  Fig.  3. 
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In  the  initial  iteration  only  the  longer  segments  of  image  1 439  with  lengths  greater 
than  30  pixels  (7.5  km)  were  considered.  Thus  in  the  test  image  (orbit  1482)  matching 
segments  of  length  n/2  =  15  pixels  and  above  were  considered  in  the  search.  The 
rather  strict  criteria  on  matching  features  were  successful  in  removing  many  erroneous 
matches,  although  it  did  also  eliminate  a  few  correct  matches  as  well.  The  subsequent 
iterations  added  progressively  fewer  matches  until  by  the  seventh  iteration  no  new 
matches  were  found.  The  primary  difficulty  lies  in  making  sure  incorrect  matches  are 
rejected  rather  than  making  sure  correct  matches  are  retained. 

The  right  panel  of  Fig.  3  shows  the  resulting  ice  displacement  vectors  for  the  image 

pa'r  aF!?‘  1 '  These  disP,acement  actors  agree  very  well  with  the  manually  generated 
and  IPAC  displacement  vectors  of  Fig.  2.  We  note  that  there  are  large  regions  where  no 
displacement  vectors  are  shown.  This  situation  arises  from  two  causes.  First,  no  ice 
displacement  vectors  will  be  generated  in  regions  where  there  are  no  lead-floe 
boundary  features  to  be  tracked.  The  second  is  that  in  regions  where  there  is  too  much 
deformation  a  boundary  segment  in  the  reference  image  will  often  be  too  badly  broken 
up  in  the  test  image  to  be  recognized.  In  the  first  case  other  features,  such  as  pressure 
ndge  patterns,  could  be  tracked  as  well  as  the  lead-floe  boundary  segments.  In  the 
second  case  the  use  of  higher  resolution  images  might  be  helpful  if  the  deformation 
were  less  severe  at  smaller  size  scales.  A  drawback  of  using  higher  resolution  pixels  in 
bAR  images  is  that  speckle  noise  becomes  a  significant  component  of  the  intensity 
variation  from  one  pixel  to  the  next. 


Fig.  3.  Search  spaces  of  long  floe-lead  boundary  segments  for  the  reference  (orbit  1439)  image  (left) 
and  test  (orbit  1481)  image  are  shown  in  the  left  and  middle  panels  respectively.  The  right  panel  shows  ice 
velocity  estimates  obtained  by  applying  the  boundary  segment  feature  tracking  algorithm  to  Fig.  1  images. 
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IV.  TRACKING  SEA  ICE  WHICH  ROTATES 


The  automated  ice  tracking  algorithms  above  work  reasonably  well  in  regions 
where  there  is  little  rotation  and  distortion,  e.g.  in  portions  of  central  pack  ice  (Fily  and 
Rothrock,  1986  and  Vesecky,  et  al.,  1988).  However  in  regions  where  there  is 
significant  ice  rotation,  >  10°  to  15°,  the  existing  image  pyramid  area  correlation 
(hierarchical  correlation)  methods  are  not  expected  to  work  well  because  the  cross 
correlation  comparisons  are  made  between  image  blocks  which  have  been  translated 
in  position,  but  not  rotated.  For  the  image  pair  of  Fig.  1 ,  upon  which  the  image  pyramid 
algorithm  works  well,  the  maximum  ice  rotation  is  some  5°.  As  noted  in  section  V,  ice 
motion  containing  shear  and  distortion  can  also  cause  errors  when  using  IPAC. 

A.  Block  Correlation  Function  in  the  Presence  of 
Rotation  and  Shear 

We  have  studied  the  effect  of  rotation  and  shear  on  block  correlation  by 
simulating  distortions  on  the  computer.  We  are  able  to  remap  an  image  by  specifying 
any  linear  transformation  between  the  coordinates  of  the  original  image  and  the  output 
image  and  then  using  bicubic  interpolation  to  resample  the  image.  In  particular,  we  can 
simulate  rotation  and  pure  shear  using  the  following  matrices: 

Rotation:  Shear: 

(cos  9  -sin  0\  /l  -tan  9  \ 

sin  0  cos  9  I  y  0  1  y 

Using  the  remapped  images  the  two-dimensional  correlation  function  for  typical  SAR 
sea  ice  scenes  can  be  calculated  under  ice  rotation  or  shear  according  to  the  matricies 
above.  For  no  rotation  or  shear,  0  =  0,  the  correlation  function  is  unity  and  no  problems 
occur.  As  9  increases,  the  correlation  decreases  and  the  probability  of  a  mismatch 
becomes  significant.  That  is,  when  we  compare  a  scene  with  itself  we  expect  unity 
correlation  coefficient;  but  we  don’t  get  it  if  rotation  or  shear  has  distorted  the  scene  we 
are  comparing  with  our  reference  scene.  One  could,  of  course,  persue  block  correlation 
with  both  translation  and  rotation.  However,  compution  time  would  increase  by  a 
factor  of  about  1 0  to  1 00. 

Vesecky  et  al.  (1987)  used  two  100x100  km  SAR  images  of  pack  ice  from  SEASAT 
orbits  1439  and  1396  to  study  the  autocorrelation  of  sea  ice  scenes  under  rotation  and 
shear.  This  experimental  study  consisted  of  extracting  twenty  non-overlapping  16x16 
subimages  at  125  meter  pixels  and  applying  both  the  above  transformations  to  each 
subimage  at  increments  of  one  degree.  The  correlation  coefficient  between  the  original 
and  the  rotated  or  sheared  image  was  then  calculated  for  a  window  size  of  1 6x1 6 
pixels,  as  used  in  the  IPAC  algorithm  of  section  II.  Fig.  4  shows  a  plot  of  the  average 
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correlation  coefficient  for  the  twenty  subimages  as  a  function  of  angle  of  rotation  or 
angle  of  shear,  0. 


5  10  15  20  25 

degrees 


Fig.  4.  Effects  of  rotation  and  shear  on  correlation  coefficient.  Averaged  two-dimensional 
cross-correlation  function  with  respect  to  angle  of  rotation  and  angle  of  shear  (9). 

If  the  correlation  coefficient  of  the  correct  match  is  low,  then  the  probability  of  a 
mismatch  increases.  In  one  case,  when  we  searched  a  local  area  around  the  correct 
match,  we  found  a  higher  correlation  coefficient  than  the  correct  one  when  the  rotation 
angle  was  1 6  degrees.  The  maximum  correlation  coefficient  was  .39  and  the  distance 
to  the  correct  match  was  about  3  pixels.  We  believe  that  rotations  of  10°  or  less  will  not 
hamper  detection  by  block  correlation. 

One  can  intuitively  interpret  the  fact  that  the  shear  curve  lies  above  the  rotation 
curved  by  considering  the  overlap  of  pixels  for  a  given  shear  or  rotation.  The  pixels  in 
the  middle  row  of  the  image  are  fixed  for  a  shear  whereas  only  the  middle  pixel  is  fixed 
for  a  rotation. 

In  summary,  the  IPAC  method  of  section  II  is  strictly  valid  only  for  translational 
motion.  However,  ice  motion  is  not  always  pure  translation,  but  has  rotational  and 
nonuniform  translational  components.  Studies  of  simulations  of  rotations  and  shears, 
Fig.  4  and  Vesecky  et  al.  (1987),  have  shown  that  IPAC  is  still  useful,  up  to  moderate 
angles  of  distortion.  In  particular,  we  believe  that  rotations  of  up  to  aboutIO0  can  be 
dealt  with  adequately  with  correlation. 

B.  Rotation  Invarient  Correlation 

Invariant  moments  are  a  pattern  recognition  technique  used  to  help  solve  the 
problem  of  identifying  objects  in  a  scene  that  changes  because  the  objects  translate  or 
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rotate  from  one  image  of  the  scene  to  the  next.  As  mentioned  above,  IPAC  will  handle 
rotations  in  sea  ice  features  up  to  about  1 0°.  Beyond  that,  rotations  will  start  to  affect  the 
ability  of  the  algorithm  to  correctly  track  features.  Rotations  of  that  magnitude  or  greater 
commonly  occur  in  places  such  as  the  marginal  ice  zone.  One  way  to  work  around 
these  large  rotations  is  to  use  a  method  that  ignores  rotations  in  a  scene  altogether. 
Invariant  moments  is  one  such  method,  whereby  invariant  moments  reduce  a  scene  to  a 
set  of  numbers,  indices  that  do  not  change  with  rotation.  We  think  that  invariant 
moments  may  also  help  sort  out  ambiguities  which  occur  with  IPAC  algorithms  applied 
to  regions  of  pack  ice  that  undergo  rotations  greater  than  the  limits  suggested  above. 

However,  the  problem  of  matching  sea  ice  features  is  more  complex  than  simply 
the  problem  of  coping  with  rotational  transformations  because  sea  ice  features  may 
change  with  translation  or  rotation  as  ice  floes  move  across  open  water  or  ice  fields. 

Floes  may  collide  with  each  other,  thereby  creating  new  features.  Rotation  of  floes  can 
affect  how  features  are  manifested  in  a  SAR  image  because  of  changes  in  aspect  angle 
(Carsey  etal.,  1983;  Curlanderet  al. ,  1987).  Translations  of  floes  might  occur  across  a 
changing  background  of  new  ice,  wind-roughened  ocean  or  pack  ice.  A  sea  ice 
tracking  algorithm  is  put  to  task  under  these  conditions. 

Although  a  variety  of  rotation-translation  invariant  measures  exist,  we  use  those 
given  by  Hu  (1962).  These  seven  indices  demonstrate  invariance  under  translation 
and  rotation.  Furthermore,  they  also  demonstrate  invariance  under  scaling  and  optical 
inversion.  Such  properties  —  invariance  under  translation,  rotation,  scaling,  and  optical 
inversion  —  make  these  indices  well  suited  for  scene  matching  problems  (Wong  and 
Hall,  1978).  To  show  why  Hu's  invariant  moments  have  these  properties,  we  will 
highlight  their  derivation. 

Hu's  invariant  moments  start  with  the  theorem  stated  by  Papoulis  (1965).  In 
essence,  this  theorem  states  that  a  scene,  over  the  (x,y)  plane  with  an  image  intensity 
function  l(x,y),  can  be  represented  by  a  unique  sequence  of  double  moments  mpq. 

These  moments  mpq  are  integrals  over  the  scene  with  integrands  of  the  form  xqyq  l(x,y). 
This  means  that  an  image  containing  many  thousands  of  pixels  large  can  be 
compressed  to  a  sequence  of  a  few  numbers,  e.g.  Hu's  invariant  moments.  Note  that 
low-order  double  moments  of  images  are  analogous  to  moments  for  objects  (mass, 
center  of  mass,  moments  of  inertia). 

Double  moments  by  themselves  show  no  invariant  properties.  However,  Hu 
derived  invariant  moments  from  double  moments  with  a  series  of  transformations.  From 
double  moments,  we  obtain  central  moments.  Under  the  translation  of  coordinates, 
central  moments  remain  invariant.  From  central  moments,  we  obtain  normalized  central 
moments.  Under  scaling,  in  addition  to  translation,  normalized  central  moments  remain 
invariant.  From  normalized  central  moments  and  with  the  theory  of  algebraic  invariants, 
we  obtain  Hu's  indices.  Under  proper  and  improper  orthogonal  rotation  and  mirroring, 
in  addition  to  translation  and  scaling,  Hu's  indices  remain  invariant.  Finally,  since  the 
magnitude  of  these  moments  may  span  many  orders,  the  typical  convention  is  to  take 
the  absolute  value  of  the  log  of  Hu's  indices. 
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Vesecky  et  al.  (1987)  have  applied  Hu's  invarient  moments  to  rotating  ice  floes  in 
the  marginal  ice  zone.  Ice  floes  which  have  rotated  90°  and  more  are  correctly 
identified  even  though  they  have  changed  appearance  somewhat  due  to  the  aspect 
angle  dependence  of  SAR  images.  However,  false  matches  were  also  encountered. 
Rotational  invarient  correlation  methods  are  subject  to  limitations  imposed  by  the 
discrete  nature  of  digital  images,  changes  in  the  background  against  which  the  floe  is 
imaged  (e.g.  surrounding  floes),  aspect  angle  dependence  of  SAR  images  and 
changes  in  the  tracked  floe  itself  (e.g.  bits  breaking  off).  In  spite  of  these  limitations  we 
think  that  rotation  invarient  correlation  will  fine  useful  application  in  tracking  sea  ice 
which  rotates.  The  essential  items  for  successful  application  of  rotation  invarient 
correlation  appear  to  be  isolating  a  tracked  floe  from  its  background  and  finding  an 
effective  similarity  measure  between  the  invarient  moments  of  two  floes. 

C.  Feature  Tracking 

The  feature  tracking  algorithm  is  inherently  robust  with  respect  to  rotation 
because  the  features  are  characterized  with  respect  to  a  coordinate  system  tied  to  the 
feature  itself  and  not  with  respect  to  geographic  or  image  coordinates.  To  understand 
this  we  consider  the  characterization  process  in  more  detail.  After  pixels  along  floe-lead 
boundaries  are  classified  and  assembled  into  segments  these  boundary  segment 
features  are  characterized  by  their  deviation  from  a  straight  line.  First,  a  'least  squares' 
criterion  straight  line  is  fit  to  a  given  boundary  segment.  The  perpendicular  distance 
from  this  straight  line  to  the  segment  is  calculated  at  regular  intervals  along  the  straight 
line.  This  set  of  deviations  characterizes  the  shape  of  a  boundary  segment. 

Recognition  of  matching  features  is  done  simply  by  comparing  the  sets  of  deviation 
measurements  for  two  candidate  features.  A  successful  match  reveals  ice  motion. 

Because  this  method  of  characterizing  a  boundary  segment  feature  is  invariant 
under  rotation  of  the  feature  one  expects  that,  in  principle,  a  feature  which  has 
undergone  rotation  will  be  as  easily  found  as  one  which  has  not.  In  the  one  pair  of  SAR 
sea  ice  images  so  far  examined  by  Vesecky  et  al.  (1988)  features  were  successfully 
tracked  using  this  algorithm  (Fig.  3).  However,  the  SAR  scenes  contained  central  pack 
ice  and  the  average  rotation  of  the  sea  ice  was  only  about  5°. 

Although  our  feature  tracking  algorithm  has  advantages  in  terms  of  following 
rotating  ice,  there  are  also  limitations  so  far  as  the  algorithm  of  section  III  is  concerned. 
The  limitation  arises,  not  in  the  characterization  of  a  feature,  but  in  choosing  what  pixels 
constitute  the  feature.  Features  are  constructed  by  first  classifying  floe-lead  boundary 
pixels  and  then  assembling  them  into  contiguous  sequences  of  pixels  having  only 
certain  allowed  classifications.  The  classification  system  uses  surrounding  pixels  for 
reference  and  thus  is  not  completely  independent  of  image  coordinates.  This  process 
allows  two  sides  of  a  lead  to  be  separate  features  and  thus  move  independently,  i.e.  it 
allows  leads  to  open  and  close.  Thus  the  lead-floe  boundary  surrounding  a  lead  may 
be  chopped  into  features  (for  tracking)  in  two  different  ways  depending  on  how  it  is 
oriented  in  the  first  and  second  images  of  an  image  pair. 
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Feature  tracking  is  a  broad  catagory  of  pattern  recognition  techniques  for 
application  to  sea  ice.  It  has  inherent  advantages  for  ice  tracking  in  the  presence  of 
rotation.  The  particular  algorithm  of  section  III  relys  on  lead-floe  boundaries  as  features 
and  thus  has  limitations  which  arises  from  the  fact  that  leads  often  open  and  close  and 
thus  the  dominant  edges  must  be  tracked  independently.  Other  sea  ice  features  do  not 
have  this  limitation.  For  example,  pressure  ridge  patterns  are  an  evident  feature  in  Fig. 

1  above.  Although  pressure  ridges  are  less  distinct  features  than  lead-floe  boundaries, 
they  are  more  stable  since  they  are  not  affected  by  opening  and  closing  of  leads.  Thus 
one  expects  that  incorporating  pressure  ridge  patterns  into  a  feature  tracking  algorithm 
would  improve  performance,  particularly  for  large  and  spatially  variable  ice  rotation. 


V.  IMAGE  RECONSTRUCTION  &  SUBTRACTION 
for  ERROR  ASSESSMENT 

A.  Image  Restoration,  Remapping  and  Image  Subtraction 

The  goal  of  this  research  is  to  produce  ice  velocity  data  for  sea  ice  science  and 
applications.  To  make  the  ice  velocity  field  estimates  of  Figs.  2  and  3  truly  useful  they 
must  be  accompanied  by  some  measure  of  the  quality  of  the  estimates.  We  have  found 
image  remapping  to  be  a  useful  technique  for  evaluating  the  quality  of  a  velocity  field. 
The  estimated  ice  velocity  field  defines  a  mapping  between  the  pair  of  images.  One  can 
use  this  mapping  to  remap  the  intensities  in  the  second  image  to  the  locations  they 
would  be  found  in  the  first  image.  It  can  be  shown  that  the  difference  between  the  first 
image  and  the  remapped  version  of  the  second  image  bounds  the  errors  in  the  motion  if 
there  is  detail  in  the  image.  Intuitively,  if  the  motion  algorithm  is  correct,  the  remapped 
intensities  subtract  perfectly  from  the  original  intensities  leaving  zero  remainder.  Where 
there  are  velocity  errors,  the  intensities  will  not  subtract  perfectly  and  the  errors  will  be 
displayed  as  bright  portions  in  the  difference  image. 

We  recognize  that  in  the  practical  world  of  SAR  remote  sensing  there  are  likely  to 
be  changes  in  intensity  between  the  two  images  in  a  pair  which  are  unrelated  to  ice 
motion.  For  example  changes  in  system  gain,  incidence  angle,  aspect  angle  and 
surface  conditions  (such  as  snow  and  wind)  can  cause  variations  in  SAR  image 
intensity.  However,  the  successful  results  of  Fig.  5  below  resulted  from  the  pair  of 
SEASAT  SAR  images  shown  in  Fig.  1.  So  we  are  optimistic  that  the  method  will  have 
application  to  a  significant  fraction  of  satellite  SAR  observations  of  sea  ice,  which  are 
less  subject  to  the  variations  mentioned  above  than  aircraft  observations. 

B.  Error  Detection  by  Image  Subtraction 

Fig.  5  shows  the  absolute  value  of  the  difference  between  the  5  Oct.  (orbit 
1 439)  image  and  the  remapped  8  Oct.  (orbit  1 482)  image  of  Fig.  1 .  The  remapped  8 
Oct.  image  was  created  using  the  smoothed  IPAC  velocity  vector  estimates  of  Fig.  2. 

The  difference  picture  has  been  displayed  so  that  image  intensity  differences  greater 
than  a  chosen  threshold  value  are  bright  and  those  below  are  dark.  Thus  high  error 
regions  are  bright  and  low  error  regions  are  dark. 
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This  image  provides  a  convenient  display  of  the  regions  where  the  ice  velocity 
estimates  are  likely  to  be  in  error.  Here  we  have  assumed  that  the  ice  itself  remains  the 
same,  as  seen  by  the  SAR,  in  the  two  images.  This  error  mao  can  be  used  in  two  ways 
First,  it  can  be  used  as  a  guide  to  the  credibility  of  the  ice  velocity  field  estimates 
provided  by  SAR.  Thus  the  ice  science  or  applications  user  is  in  a  position  to  know 
when  and  where  the  ice  velocity  data  provided  to  him  are  dependable  or  likely  to  be  in 
Qffor.  Second,  the  error  data  can  be  used  as  feedback  within  an  ice  velocity  estimation 
algorithm  to  improve  the  velocity  estimates. 


Fig.  5.  Image  reconstruction  and  subtraction  map  based  on  the  ice  velocity  field  of  Fig  2  Regions 
where  the  ice  velocity  estimates  are  likely  to  be  in  error  are  bright  and  low  error  regions  are  dark. 


C.  Interpretation  of  Image  Subtraction  Results 

By  comparing  the  error  display  of  Fig.  5  with  the  ice  velocity  field  as  shown  in 
Fig.  2  we  find  that  we  can  make  an  interesting  interpretation  of  the  error  display.  First, 
we  note  that  the  large  dark  areas  where  errors  are  low  correspond  to  large  blocks  of  ice 
which  move  as  rigid  plates.  These  large  units  were  also  pointed  out  by  Fily  and 
Rothrock  (1 986).  Second,  we  see  that  the  large  error  (bright)  regions  tend  to  fall  along 
lines  which  separate  the  large  rigid  plates.  It  is  in  these  regions  that  the  large  plates 
aoparently  grind  against  one  another  producing  velocity  shear  which  deforms  the  ice. 
The  deformation  of  the  ice  in  these  regions,  coupled  with  the  fact  that  the  ice  velocity 
estimates  are  discrete  and  smoothed,  causes  errors  in  ice  velocity  estimates.  Thus  we 
suggest  that  the  remapping  error  plots  such  as  that  shown  in  Fig.  5  can  be  used  to 
segment  an  ice  velocity  field  into  dark  regions  which  tend  to  be  large  rigidly  rotating 
plates  and  bright  regions  where  there  is  much  shear  and  ice  deformation. 


61 


VI.  SYNTHESIS  of  a  GENERAL  SEA  ICE  TRACKING  ALGORITHM 


In  work  described  here  as  well  as  that  by  Fily  and  Rothrock  (1986),  Michael  Collins 
(1 987),  Tom  Hiroshi  and  others  we  find  several  schemes  for  automated  sea  ice  tracking 
and  the  estimation  of  sea  ice  velocity  fields.  At  this  stage  it  appears  that  the  different 
algorithms  each  have  their  own  relative  advantages  and  disadvantages.  The  relative 
advantages  of  IPAC  and  feature  tracking  have  been  discussed  by  Vesecky  et  al.  (1988). 
No  one  algorithm  is  clearly  superior  over  a  wide  range  of  sea  ice  conditions.  Thus  it 
appears  that  a  synthesis  of  several  methods  will  be  most  likely  to  provide  a  robust, 
automated  sea  ice  tracking  algorithm. 

To  illustrate  the  roles  of  various  methods  we  briefly  describe  several  techniques 
and  where  they  might  fit  in  an  algorithm  synthesis.  Collins  (1987)  has  shown  that  a 
matched  filter  approach  implies  that  the  image  data  should  be  filtered  before  any 
correlation  algorithms  are  applied.  The  resulting  ’high  pass'  filter  would  emphasize  the 
edges  in  the  images  and  make  correlation  work  more  effectively.  Before  IPAC  or 
related  block  correlation  techniques  can  be  applied  the  amount  of  ice  rotation  shoule  be 
assessed  and  compensated  for  if  greater  than  about  10°.  Feature  tracking  and  rotation 
invarient  correlation  could  be  used  to  estimate  ice  rotation  (so  that  it  could  be 
compensated  for  in  an  IPAC  scheme)  or  provide  alternatives  to  block  correlation 
techniques.  The  point  here  is  that  the  different  ice  tracking  methods  can  be  used  to 
mitigate  the  deficiencies  in  one  another. 

On  way  of  thinking  about  the  structure  of  a  synthesis  algorithm  is  to  consider  series 
and  parallel  strategies.  In  a  series  strategy  the  various  component  algorithms  would  be 
applied  in  sequence.  An  example  of  this  type  of  algorithm  was  presented  by  Curlander 
(1987).  The  idea  is  to  use  a  feature  tracking  algorithm  first  to  define  the  gross  ice 
motion,  including  rotation,  if  any.  This  preliminary  estimate  of  the  velocity  field  would 
then  be  used  as  a  basis  for  application  of  an  IPAC  type  technique.  At  each  location 
where  block  correlation  is  to  be  carried  out  the  scene  in  the  test  image  (second  in  time) 
would  be  counterrotated  according  to  the  rotation  estimate  of  the  first  step  so  that  block 
correlation  would  effective  as  discussed  in  section  IV  above. 

An  example  of  a  parallel  structure  is  shown  in  Fig.  6  (Vesecky,  1987).  Here  IPAC  is 
the  primary  technique  with  feature  tracking  and  rotation  invarient  correlation  being 
employed  to  check  consistency  and  in  locations  where  IPAC  yields  low  confidence 
estimates.  Referring  to  Fig.  6  we  see  that  environmental  data,  primarily  surface  winds, 
are  used  to  select  the  pair  of  SAR  scenes  which  are  the  input  to  the  tracking  algorithms. 
Feature  tracking  is  used  as  subordinate  technique  to  confirm  or  disagree  with  the  IPAC 
estimates.  At  each  level  in  the  pyramid  the  ice  velocity  field  estimated  by  IPAC  is 
checked  for  confidence  (as  discussed  in  section  II. B)  and  consistency  with  the  feature 
tracking  (discussed  in  section  III)  velocity  estimates.  Rotation  invarient  correlation 
methods  could  come  into  play  here  to  revise  IPAC  velocity  estimates  which  are  either 
low  confidence  or  inconsistent  with  feature  tracking  estimates.  Velocity  estimates  which 
remain  below  a  confidence  and  consistancy  threshold  are  then  removed  and  a 
substitute  generated  by  interpolation  from  surrounding  velocity  vectors  (see  section 
II. C).  After  the  lowest  level  in  the  IPAC  pyramid  has  been  completed  error  estimates  are 
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computed  by  using  the  image  reconstruction-subtraction  scheme  (section  V).  After  an 
unspecified  last  effort  to  resolve  vectors  which  are  likely  to  be  in  error  the  system 
produces  its  output,  namely:  V(x.y),  sea  ice  velocity  estimates  as  a  function  of  time  and 
geographic  position  (Figs.  2  and  3)  and  error  estimates  for  V(x,y)  as  shown  in  Fig.  5. 


COMPREHENSIVE  SAR  ICE  TRACKER 


SAR  I  mag*  Data 


Vfc#(i.y)  V  Error  Estimai* 


Fig.  6.  Schematic  block  diagram  for  a  comprehensive  automated  system  for  estimating  sea  ice  velocity 
from  SAR  images.  This  scheme  incorporates  several  techniques  by  using  a  parallel  algorithm  structure. 


VII.  DISPLAY  OF  SEA  ICE  VELOCITY  DATA 

While  accurate  ice  velocity  estimates  are  most  important,  effective  display  of  the 
results  to  the  user  is  a  vital  link  in  using  SAR  images  in  sea  ice  science  and 
applications.  For  it  is  the  user  who  will  interpret  the  ice  velocity  fields.  An  ineffective 
display  scheme  can  not  only  foil  creative  interpretation  by  producing  artifacts  in  the 
display,  but  can  also  obscure  significant  features. 

As  an  example,  consider  the  commonly  used  line  or  arrow  display  of  Figs.  2  and  3. 
While  this  display  seems  perfectly  adequate,  note  that  artifacts  which  draw  the  eye  are 
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produced  when  the  lines  happen  to  be  horizontal,  vertical  or  inclined  at  45°,  i.e.  when 
the  arrow  heads  follow  the  tails  of  their  neighbors.  Clearly  this  artifact  is  related  to  the 
choice  of  coordinates  and  has  no  physical  significance  in  ice  science. 
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Fig.  7.  Hexagon  display  scheme  for  emphasizing  shear  and  distortion  in  a  sea  ice  velocity  field. 


Other  types  of  displays  are  possible  and  we  think  can  contribute  significantly  to  the 
interpretation  of  sea  ice  velocity  data.  For  example,  suppose  that  one  wants  to 
emphasize  the  distortion  or  shear  in  an  ice  velocity  field.  The  scheme  illustrated  in  Fig. 

7  uses  a  pattern  of  octagons  placed  on  the  original  SAR  image  at  regular  intervals.  The 
ice  velocity  field  is  then  allowed  to  move  each  vertex  of  each  octagon.  If  the  ice  velocity 
field  is  uniform,  an  octagon  retains  is  shape.  If  there  is  velocity  shear  or  other  distortion, 
then  each  octagon  distorts  accordingly.  The  velocity  field  used  to  construct  Fig.  7  is  that 
of  Fig.  2.  We  have  darkened  the  hexagons  which  are  distorted  to  emphasize  regions  of 
shear  and  distortion.  Note  the  correspondence  between  Figs.  5  and  7  which  both  make 
clear  the  large,  rigidly  moving  plates  of  ice  and  the  regions  of  shear  and  distortion 
between  them,  where  they  grind  together. 


VIII.  CONCLUSIONS 

1 .  The  image  pyramid  area  correlation  algorithm  (IPAC)  provides  a  useful,  though 
not  ideal,  means  of  automating  the  process  of  deriving  sea  ice  velocity  information  from 
SAR  images. 

2.  Feature  tracking  using  floe-lead  boundary  segments  produces  sea  ice  velocity 
estimates  which  agree  well  with  IPAC  and  manual  results,  but  are  limited  to  regions 
where  features  are  found.  Extension  to  pressure  ridges  and  other  features  is  indicated. 

3.  Current  block  correlation  schemes  are  likely  to  be  in  error  when  ice  rotation  is  > 
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10°  between  images.  Feature  tracking  and  rotational  invarient  correlation  can  provide 
tracking  capability  when  the  ice  rotates. 

4.  Image  reconstruction  and  subtraction  provides  a  convenient  measure  of  regions 
where  ice  velocity  estimates  are  likely  to  be  in  error. 

5.  A  robust,  automated  sea  ice  tracking  algorithm  will  need  a  variety  of  techniques 
to  operate  successfully  over  a  wide  range  of  sea  ice  conditions.  A  synthesis  of  several 
available  techniques  is  indicated. 

6.  Effective  display  of  sea  ice  velocity  information  is  important  for  effective 
interpretation  of  the  data. 
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remote  sensing  of  sea  ice  motion  using  floe  edge  and  pressure  ridge  features  in 

SAR  IMAGES 
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Synthetic  aperture  radar  (SAR)  provides  an  excellent  means  of  observing  the 
movement  and  distortion  of  sea  ice  over  large  temporal  and  spatial  scales.  SAR 
observations  are  not  affected  by  clouds  or  lack  of  sunlight.  The  European  Space  Agency’s 
ERS-1  satellite  as  well  as  others  planned  to  follow  will  carry  SAR*s  over  the  polar  regions 
beginning  in  the  early  1990's.  Ground  stations,  such  as  NASA's  Alaskan  SAR  Facility,  will 
interpret  these  SAR  images  in  terms  of  geophysical  measurements  of  sea  ice.  A  key 
component  in  utilization  of  SAR  image  data  is  an  automated  scheme  for  extracting  the  sea 
ice  velocity  field  from  a  sequence  of  SAR  images  of  the  same  geographical  region.  The 
process  of  retrieving  motion  information  from  a  sequence  of  SAR  images  revolves  around 
two  principal  items,  namely: 

1.  Identifying  regions  in  the  image  of  a  sea  ice  field  to  use  as  lie  points'  for  ice 

floe  tracking. 

2.  Finding  the  same  tie  points  in  subsequent  images  of  the  sequence. 

Once  the  ice  displacement  is  known  the  ice  velocity  can  easily  be  calculated  since  the 
elapsed  time  between  the  SAR  images  is  known.  Two  general  techniques  for  automated 
sea  ice  tracking  have  emerged  from  recent  research,  namely:  image  pyramid  area 
correlation  (hierarchical  correlation)  and  feature  tracking.  The  first  method  uses  a  pyramid 
(variable  resolution)  data  structure  and  relies  on  finding  the  maximum  block  correlation 
between  areas  in  the  two  images.  It  performs  in  a  computationally  efficient  manner  when 
the  ice  being  tracked  rotates  less  than  10°  to  15°.  For  ice  which  rotates  more  than  15° 
between  SAR  observations,  such  as  in  the  marginal  ice  zone,  the  feature  tracking  method 
appears  preferable.  Whilst  image  pyramid  techniques  have  been  reported  elsewhere,  we 
focus  here  on  feature  tracking.  Distinctive  features  are  located  and  characterized  in  each 
image;  thus  converting  each  image  from  a  set  of  pixels  to  a  set  of  features.  These  features 
are  then  used  as  tie  points  to  estimate  ice  motion.  Since  each  image  is  converted  from 
millions  of  pixels  to  sets  of  tens  of  features  each  characterized  by  about  a  hundred 
numbers,  the  computation  time  required  to  establish  matching  features  is  tractable.  To 
observe  ice  motion  when  the  ice  rotates,  the  features  must  be  characterized  in  a  manner 
which  is  independent  of  orientation.  Here  we  consider  two  types  of  features:  floe-lead 
boundaries  and  pressure  ridges.  In  both  cases  we  process  the  SAR  image  data  to 
emphasize  the  feature  under  consideration.  For  floe-lead  boundaries  the  image  is 
segmented  into  ice  and  water  areas  by  thresholding  and  boundary  segments  identified  by 
simply  looking  for  ice  pixels  next  to  water  pixels.  Pressure  ridges  are  emphasized  by 


PmctetBnp  of  1GARSS  88  Symposium,  Edinburgh,  Scotland,  13-16  Sept.  19S&  Ref.  ESA  Sf-284  (IEEE  88CH24974). 
Published  by  ESA  Publications  Division.  August  1988. 


* 


F 

418 


i  F  VESECKY  SlAL 


expanding  the  brighter  portions  of  the  image  and  thresholding.  Gaps  in  filamentary 
features  are  filled  in  using  a  growth  and  erosion  technique.  Boundary  segments  are 
characterized  by  fitting  a  straight  line  to  the  segment  and  noting  the  deviation  of  the 
segment  from  this  line.  Pressure  ridges  are  characterized  by  using  a  fitted  line  as  the  trunk’ 
of  a  tree  structure  with  branches  jutting  out  from  the  trunk  at  different  points  and  angles  for 
different  distances.  Since  these  features  are  referenced  to  a  line  and  the  line  can  have  any 
location  and  orientation,  the  set  of  numbers  characterizing  the  feature  is  independent  of  the 
feature’s  location  and  orientation  in  the  image.  Once  each  image  is  converted  into  a  set  of 
features,  the  features  are  matched  by  correlation.  The  iterative  algorithm,  used  in  the 
search,  looks  for  high  confidence  matches  between  the  most  distinctive  features  and  uses 
these  to  guide  further  searches  in  terms  of  location.  We  apply  these  methods  to  several 
pairs  of  Seasat  SAR  images  of  sea  ice  from  the  Beaufort  Sea  and  present  the  results  The 
results  compare  favorably  with  ice  motion  estimates  made  by  visually  finding  tie  points  in 
the  same  images  using  a  split  screen  video  display. 


